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Abstract 

There are always differences between theoretical and experimental results in the study of lifting 
surfaces. Bounding box control volume measurements infrequently yield exact conservation of mass 
or consistent values for lift and drag coefficients. Numerically calculated wakes often differ from 
experimental data. Quite often, an empirical correction can be applied to fit theory to experiment 
to account for these differences. However, as the demands for state of the art foil design increase, 
fluid dynamicists are pressed to look carefully at these inconsistencies in order to improve current 
design and analysis methods. Using a Reynolds Averaged Navier Stokes (RANS) computer code and 
a highly refined fluid mesh, one can begin to explore the subtle characteristics of the fluid flow in 
the entire domain and the details of certain key regions around a foil. Specific areas of great interest 
are: flow around the trailing edge, flow within the boundary layer, wake profiles and the influence 
of tunnel wall boundaries in experimental facilities. 

The overall goal of this thesis is to resolve some of the discrepancies between theoretical results 
and experimental data. A computer code has been developed to generate the geometry for the fluid 
flow domain surrounding an arbitrary foil shape at a specified angle of attack in the MIT Marine 
Hydrodynamics Laboratory (MHL) water tunnel. This geometry is provided as input data for the 
RANS solver. A suite of software tools are developed to provide post processing analysis to compare 
the RANS solution with other numerical techniques and experimental measurements. 

Through the use of case studies, the numerical results of the RANS code are compared with 
recent MHL experimental data and other computational tools. A comparison is made between 
the experimental and RANS code results using a control volume analysis. Boundary layer and 
wake profiles are also compared. A correction scheme is developed to extrapolate experimental 
measurements to unbounded fluid flow. 
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Mathematic 


: Notation(continued) 


5 


domain of foil surface 


t 


tangent unit vector 


tx 


J-comp. of tangent vector on 


ty 


y-comp. of tangent vector on 


T 


fluid shear stress 


U inf 


Freestream velocity 


u 


local ^-velocity 


\Jl~ns/p 


Friction Velocity 


V 


local y-velocity 


V 


velocity 


r 


dimensional circulation 


7 


vorticity 


76 


bound vorticity 


7/ 


free vorticity 


P 


dynamic viscosity 


V 


kinematic viscosity 


P 


density of fluid 


V 2 


Laplace Operator 



surface 

surface 
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Chapter 1 

Introduction 



1.1 Overview 

Recent application of advanced computational methods to design new hydrofoil 1 sec- 
tions for propeller blades has yielded inconsistent results. There is a need in the hy- 
drodynamic research community to undertake an in depth study of two-dimensional 
lifting surfaces since these sections are the fundamental building blocks of most three 
dimensional design methods. Additionally, a sound methodology should be applied 
to make effective comparisons between numerical and experimental results. Once the 
differences between numerical predictions, experimental measurements and real appli- 
cation are understood, those lessons can be incorporated into computational methods. 
Some improvement will come from better understanding of fluid flow within boundary 
layers near the foil surface. Other areas in need of refinement are the effects of wake 
diffusion and flow separation near the trailing edge. 

Ultimately, yielding these improvements will be a multi-step process that involves 
significant computational and experimental effort far beyond the scope of a single 
thesis. The work herein represents the first few steps towards the goal of developing 
a more complete understanding of the subtleties of fluid flow around two-dimensional 
foils. 

'Throughout this thesis, the words hydrofoil, 2-D lifting surface and foil will be used as synonymous terms. 
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1.2 Background 

1.2.1 Why Is There Interest in 2-D Flow? 



Presently, the hydrodynamic research community has turned its attention back to 
the study of two-dimensional lifting surfaces. The renewed interest is derived from 
recent attempts to design advanced hydrofoil sections for propeller applications that 
demonstrate improved cavitation characteristics without incurring a significant drag 
penalty. 

Recent work in the MIT MHL[17] water tunnel by Jorde[16] and Kimball[19] on 
foils with non-traditional camber distributions and unique trailing edges indicate 
that the above goals may be attainable. However, they also point out that our 
knowledge of two-dimensional flow around foils is incomplete. Bloch [7] developed 
a methodology for adding corrections to inviscid propeller design codes to account 
for anti-singing trailing edges. Several issues have come to light. First, the current 
computer design codes work well for conventional foils. However, the current codes do 
not work well for foils with cupped or blunt trailing edges or for foils with advanced 
camber distributions designed to delay cavitation inception. 

1.2.2 A Recent Design Problem 

As part of an advanced technology demonstration(ATD), new high speed propellers 
were designed and tested at the Carderock Division, Naval Surface Warfare Center. 
The design effort required to achieve the improved propeller performance using ad- 
vanced blade sections was monumental. The design process involved iterating several 
times between using propeller design computer codes and model testing[3]. The costs 
and time associated with this type of design procedure are prohibitive. 

It is not practical to expect that any commercial operator buying a propeller 
for a ship could fund such a research and development effort. But, the commercial 
operators do want the advantages these new foil sections offer. The way to make 
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designing with advanced foil sections cost effective is to eliminate the costly iteration 
between design computer codes and model testing by developing a database that 
generalizes these new geometries as suggested by Bloch. Prudent application of model 
testing and RANS analysis for a broad geometry of foil shapes could provide the 
necessary data to formulate correction routines for standard propeller design codes 
without significantly increasing the computation time for a converged solution. 

1.3 Motivation for Modeling 2-D Foils With RANS 

An incompressible Reynolds Averaged Navier Stokes (RANS) solver was chosen for 
the bulk of the numerical analysis in this thesis. A RANS computer code solves 
the equations of motion for each fluid element throughout an entire domain. When 
compared to other solution methods, such as an inviscid panel method, RANS is very 
time consuming and requires large amounts of random access memory (RAM) and 
central processing unit (CPU) time. This is a significant disadvantage, especially 
when one wants to analyze several foils at several angles of attack and Reynolds 
numbers. But, there are many aspects of a RANS solution that are attractive. 

The voluminous output from a RANS computer code contains all of the flow char- 
acteristics for every fluid element throughout the entire domain such as: Velocity, 
Pressure, and Shear Stress. Provided the flow domain is discretized with sufficient 
resolution, you can also extract and measure subtleties of the flow like boundary layer 
velocity profiles and thickness, and flow separation under adverse pressure gradients. 
It is desirable to be able to study the characteristics of a viscous wake behind a hy- 
drofoil. Also, with the abundance of data available in the solution, it is easy to make 
a direct comparison between conditions calculated at a prescribed location in the flow 
field and those measured in a geometrically similar experiment. 

A large part of this study involves characterizing the differences between numerical 
solutions and experimental measurements. In an unbounded flow, such as a 2-D 
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foil fixed in a uniform stream, RANS codes provide a good validation for computer 
programs such as PAN2D [21] and XFOIL [10] which are both inviscid panel methods 
coupled with integral boundary layer solvers. Neither PAN2D nor XFOIL is currently 
capable of modeling a foil in a flow constrained by walls including the viscous effects 
of the boundary layers, which is the case in water tunnel experiments. A RANS 
code can serve as a liaison between unbounded numerical codes and the bounded 
case of experiments. For an equivalent foil geometry, the RANS code can be used to 
characterize the differences in lift, drag and other properties for both bounded and 
unbounded flows. A methodology for doing this is presented in Chapter 5. 

1.4 Objectives 

As was alluded to in Section 1.3, one overall goal in this thesis is to provide a scheme 
for feeding back results and measurements taken in experiments as corrections that 
can be applied to the computationally efficient inviscid panel methods. Several steps 
need to occur before this goal can be realized in a substantial way. Methods need 
to be developed to efficiently generate input files for the RANS solver. As will be 
shown in Chapter 3, accurate geometric representation of a flow domain is perhaps 
the most demanding part of obtaining the solution. Accordingly, a great deal of effort 
was devoted to ensuring that the RANS domain geometry was exactly the same as 
the inviscid computer codes and the experimental setups. Computer codes need to 
be developed to reduce the output data to a tractable and meaningful form. The 
following are the critical path or the enabling objectives used to achieve the above 
goals: 

1. Review recent MIT MHL Water Tunnel experimental results for advanced section 
two-dimensional hydrofoils. 

2. Write a computer code that generates the fluid flow geometry input files for the 
RANS solver. 
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3. Analyze foils in unbounded flows using the RANS code and PAN2D. 

4. Analyze foils in bounded flows using the equivalent geometry of the MIT MHL. 

5. Develop corrective factors to apply to the inviscid solvers which are based on the 
differences found between the RANS and experimental results. 

6. Demonstrate how RANS solutions can be used to formulate experimental test plans. 

1.4.1 Rapid Generation of Flow Geometry 

To support ongoing research in the MIT Marine Hydrodynamics Laboratory, a two- 
dimensional lifting surface analysis tool is developed in this thesis to accurately model 
the fluid flow around hydrofoils in the water tunnel. It is very useful to have good 
predictions of a hydrofoil’s performance characteristics prior to conducting an exper- 
iment. Thorough empirical evaluation of hydrofoils requires taking many measure- 
ments with varying geometry and Reynolds numbers. It is common to take measure- 
ments for a single foil geometry at several angles of attack and Reynolds number. 
Therefore, one requirement for any computational tool used in the MHL is that it 
be easy to change the foil angle of attack and other test conditions such as scale and 
Reynolds number. Part of the work in this thesis is the development of the computer 
code FIT2D(Foil In Tunnel, Two-Dimensional). FIT2D is an interactive fluid mesh 
geometry generator. The user can arbitrarily specify: angle of attack, grid resolu- 
tion, and Reynolds number as well as many other lesser parameters. The output files 
generated are used as input for a Poisson equation (V 2 0 = Const) grid refinement 
program. After the grid is refined, it becomes input for a RANS solver for analysis. 

1.4.2 Resolving Differences Between Experiments and Numerics 

In sections 1.3 and 1.2.1, it is asserted that to make the current foil design computer 
codes work better, corrections should be incorporated to account for the physical 
effects that are not modeled in computationally efficient inviscid solutions. Kimball 
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found that there is uneven boundary layer growth on tunnel walls due to the presence 
of a foil that is generating lift [19] . This difference between the upper and lower walls 
affects the flow and the resulting lift and drag measurements. So there are really 
two sides to the correction scheme. One is using RANS to quantify the effects that 
the tunnel walls have on the foil lift and drag measurements as a result of uneven 
boundary layer growth and the potential flow imaging effect. The other is using the 
experimental measurements to apply corrections to computer codes for physics that 
are not captured by the numerics such as re-attachment of separated flows, transition 
back to a laminar boundary layer, and vortex shedding phenomena. 

1.5 Organization 

In Chapter Two, all of the relevant theory of two-dimensional lifting surfaces is out- 
lined. It provides the mathematical statement of the lifting problem. Similarly, 
Chapter Three is a detailed description of the numerical aspects of the computa- 
tional solution. 

Chapter Four outlines the methodology for post processing all of the RANS data. 
It provides an overview of the suite of programs that were developed and used to 
calculate lift and drag, and analyze wake profiles and boundary layers. 

Chapter Five compares and contrasts the RANS output to other numerical meth- 
ods and experimental measurements. Results of case studies for two different foils are 
presented. Figures, Graphs, and Tables are used to demonstrate the differences and 
similarities that occur. A methodology is proposed for applying corrections hydrofoil 
computer programs to gain better agreement with precise experimental results and 
real applications. A demonstration of how RANS solutions can be used to formulate 
experimental test plans is presented. 

Chapter Six is a summary of results and conclusions. Future topics and directions 
for research in this area are discussed. 
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Chapter 2 

Theory 



2.1 Overview 

This study focuses on 2-D fluid flow, therefore it is appropriate that this chapter 
outlines and develops the pertinent fluid mechanics theories for lifting surfaces and 
boundary layers. These theories have been previously derived in many texts and 
technical papers. The developments and derivations contained in this thesis represent 
the specific application to the problem of thin non-cavitating hydrofoils in bounded 
and unbounded incompressible two-dimensional flow domains. 

First, the basic lifting problem is defined as a boundary value problem and the 
geometry is specified. Subsequent sections draw attention to aspects of lifting surface 
theory that are fundamental to the computational methods employed and to the 
physical phenomena that are modeled. 

2.2 Two-Dimensional Lifting Flows 

Lifting surfaces have many applications in marine hydrodynamics. Two dimensional 
foil sections are used extensively in the design of rudders, keels, skegs, propellers, and 
other appendages, such as ducts, where a desired force is generated normal to the 
onset flow. 

This study is limited to a “typical” lifting surface as described in Newman [24] 
where the thickness is much smaller than the chord length (t(x)<g.C). Further, only 
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incompressible two-dimensional flows are considered where the span of the foil is of 
infinite length and the resultant flow over the surface is chord-wise in direction. Figure 
2-1 shows how the foil geometry is defined with respect to the coordinate axes and 
the onset fluid flow. The foil is fixed in a stationary reference frame. The velocity 




Figure 2-1: 2-D Foil Coordinate System 

field ( Ui n f ) propagates in the direction of the positive x-axis, and accordingly the 
leading edge of the foil points opposite to the onset flow. Although the pivot point is 
arbitrary for an unbounded flow, it becomes important when walls are present close 
to the foil surface. In this work, care was taken to specify the exact x,y location of 
the pivot point for comparing numerical results to those obtained in the MHL water 
tunnel. The angle of attack (a) is defined as the angle that an imaginary line drawn 
from the leading edge to the trailing edge makes with the x-axis. Here it is assumed 
that the “nose to tail” line intersects the fore and aft end of the mean camber line 
at the zero angle of attack. If the trailing edge is blunt, beveled or has some other 
treatment such as a splitter plate, the mean camber line ends at the midpoint of a 
line drawn from the upper and lower points of the trailing edge surface. The mean 
camber line is defined as: 

y m {x) = |(y«(ar) + yi(x)) (2.1) 
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Foil thickness is defined as: 



Vt{x) = y u (x) - yi{x) (2.2) 

Steady state flow is assumed for all cases. The Kutta Condition[l8] applies at the 
trailing edge requiring that the fluid velocity be finite. With the above limitations 
and the stated geometry, the following boundary value problem results: 



VV = 0 (2.3) 

V = 0,on y u (x) and yi(x) (2.4) 

V</> < oo, at the trailing edge, Kutta Condition (2.5) 

V = 0,on walls if present (2-6) 

Vcf) = Uin/ i, as x, y — >• oo (2.7) 



Equations 2.3 - 2.7 are the complete mathematical statement of the lifting problem for 
a stationary foil with onset flow. Two types of solutions to this problem are pursued. 
A Reynolds Averaged Navier Stokes solver is used to solve the steady state viscous 
flow around the foil. The other method is a potential flow vortex lattice method 
where the no-slip rigid boundary conditions of Equations 2.4 and 2.6 are relaxed. In 
inviscid potential flow, the wall and foil surface boundary conditions become: 

V • n = 0 (2.8) 

Now that the boundary value problem is defined for the viscous and inviscid problems, 
the solutions can be formulated. 

2.3 Linearized Two-Dimensional Theory 

The use of linear potential theory is very powerful for analyzing lifting flows. A simple 
example of a lifting potential flow is a point vortex fixed in a uniform stream as shown 
in Figure 2-2. The linear potential for a free stream with a point vortex located at 
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3'0't Vo IS. 

$ = U inf x + Im^-(log(x -x 0 + i(y - y 0 )) (2.9) 

This potential is a solution to the Laplace equation (2.3) throughout the flow field 
except at the location of the point vortex. Using the principle of superposition, a 
more complicated flow can be formulated by adding a group of simple flows together 
that satisfy the kinematic boundary conditions and Laplace. 




Figure 2-2: Point Vortex in Uniform Stream 



2.3.1 Lift 

In potential flow the lifting force acting on a foil is obtained using the Kutta-J oukowsky 
theorem, which states that for any two-dimensional body, moving with constant ve- 
locity in an unbounded inviscid fluid, the hydrodynamic pressure force is directed 
normal to the velocity vector and is equal to the product of the fluid density, body 
velocity and the circulation about the body[24]. The mathematical statement of this 
theorem applied to thin foils is: 

L = (p — Poo)dx = pUju dx = pUT (2.10) 

In a potential flow that is formulated by a distribution of discrete point vortices, 
application of 2.10 is a simple matter of summing up the strengths of the point 
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vortices enclosed by a simple contour. 



2.3.2 Drag 

In a pure inviscid potential flow solution, there is no pressure or form drag. However, 
when fluid viscosity is added to the solution, a viscous boundary layer(section 2.7) is 
formed near the foil surface. Drag on the foil is present as a result of two phenomena. 
First, shear stresses in the boundary layer cause skin friction drag. Second, due to the 
presence of the boundary layer, the pressure recovery at the trailing edge is incomplete 
which causes a drag force. The pressure drag is hard to calculate or measure because 
the pressure changes that cause the drag are isolated to a small area near the trailing 
edge of the foil. Also when compared to the lift, drag is usually a very small quantity. 
The errors associated with computing the drag are nearly equal to the drag itself. 

Two other methods of computing the drag are presented later in sections 2.5.1 and 
2.5.2. These methods are not as accurate as direct integration of pressure and shear 
stress but provide additional verification of the integrated quantities. 

2.3.3 Lift and Drag by Pressure Integration 

When the viscous effects are included, Kutta-Joukowsky does not apply for calculating 
lift. The drag is no longer zero. However, if the pressure and shear stress in the fluid 
adjacent to the surface(S) are known, the lift and drag can be obtained by direct 
pressure integration on the foil surface around the perimeter of the foil. The lift and 
drag forces on the foil due to pressure integration are: 

F Dp = J s P-n x dS (2.11) 

and 

F L = J^P -n y dS (2.12) 

where S is the foil surface and P is the fluid pressure at the foil surface. In addition 
to the normal pressure to the surface, there is a viscous shearing stress (r„ s ) acting 
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tangent to the surface. In a similar fashion to the above equations, shearing forces can 
be calculated by integrating around the foil and resolving into x and y components. 
Typically, the y component of the shearing forces is neglected because the integral 
quantity of these forces is insignificant compared to the y component of the pressure 
forces. The drag force on a foil due to viscous shearing forces is: 

F Dv — J^Tns ■ t x dS (2-13) 

The total drag on the foil is the sum of the pressure and viscous drag terms, or: 

Fd = Fop + F^y (2-14) 

2.3.4 Viscous Effects on Lift 

A detailed overview of the viscous effects on lift is contained in Coney [9]. Some 
pertinent aspects of his discussion are presented here to provide needed clarity. As 
shown in Figure 2-3 1 , a foil in the presence of a real fluid flow will have a thin viscous 
boundary layer formed on its surface. This, in essence, alters the effective thickness 
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Figure 2-3: Viscous Boundary Layer on 2-D Foil 

and camber distribution of the foil. The camber is reduced on the aft portion of 
the foil. In general, for the same foil geometry, the difference between the potential 
flow and viscous flow solution is that the change in camber distribution will cause an 
apparent reduction in lift coefficient. Boundary layer characteristics are covered in 
more detail in section 2.7 

'Figure 2-3 created by Scott D. Black 
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2.4 Vortex Lattice Methods(VLM) 



The basic concepts of vortex lattice theory are contained in Newman [24] and Ker- 
vvin [18]. In light of the long computation time to obtain a RANS solution, it is a 
prudent step to incorporate a vortex lattice analysis code into the computer code that 
generates the geometry for the RANS solver. This serves two purposes. First, it pro- 
vides an initial estimate for the lift coefficient(C)) and the circulation distribution^) 
on a foil. Second, it provides the dividing streamline in the wake by extraction of 
the field point velocity due to the circulation and free stream potential. This is de- 
sirable because the numerics of the RANS solver works better if the fluid domain 
discretization is adapted to the location and trajectory of the wake.(See Section 3.7) 

Vortex lattice methods are simple and efficient. The earliest vortex lattice method 
is attributed to Falkner[12] as referencedby Kerwin[18]. The derivation presented here 
is analogous to Kerwin’s method. The primary difference is that the actual mean line 
of the foil is used vice the linearized foil surface. 

The formulation of the vortex lattice method is basic. The mean camber line of 
a given foil is determined using equation 2.1. The camber line is broken up into 
a discrete number of elements of panels. The panels are spaced using a “cosine” 
method where the panels at the leading and trailing edge of the foil are smaller than 
the middle panels. A point vortex of unknown strength is located at the mid-point 
of each of the m panels. Control points where the kinematic boundary conditions are 
imposed are located at the end of each panel. Figure 2-4 is an example of a vortex 
lattice distributed on a mean camber line. 

2.4.1 Cosine Spacing 

Cosine spacing is commonly used in many lifting surface and vortex lattice methods 
as the spacing algorithm for distributing point vortices and control points. It is a 



28 



0.6 
0.5 
0.4 
0.3 

0.2 

■s 

2 . 01 1 - 
o 
> oF- 

-0.1 r 

-0.2 r 
-0.3 r 
-0.4 r 



u 



inf 



FP P B — B B B B — B 



Point Vortex 
Control Point 



J— I I UlJ L_ 



0.25 



0.5 

x/Chord 



-I I I I 1 I 1 !_ 



0.75 



Figure 2-4: Vortex Lattice on a 2-D Mean Line 



simple method where an auxiliary variable s is defined such that: 

Slen , 



S = 



•(1 — cos s) 



(2.15) 



where s = 0 is the leading edge of the mean line and s = 7T is the trailing edge. In 
the auxiliary coordinate system, the interval of the arc is divided into N equal panels 
each of s = n/N length. To establish the locations for the point vortices and the 
control points the following equations are used: 



5 4 n ) = -TpU - cos ( 



( n — l/2)7r 

N 



)] 



s c( n ) = %4l -cos(^)] 

2.4.2 Obtaining the Vortex Lattice Solution 



(2.16) 

(2.17) 



The boundary condition on the foil is prescribed by equation 2.8. The strength of 
the point vortices are to be determined such that this boundary condition is satisfied 
at each of the control points. First it is instructive to have the expression for the 
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velocity field induced at a point x,y by a point vortex located at the origin[24]. 



V = 



T iy-jx 



(2.18) 



2 7T (a; 2 + y 2 ) 

This equation can easily be converted to the velocity field induced by a point vortex 
located at x v ,y v on a control point located at x c ,y c by making a simple substitution 
of the difference in the coordinates for x and y. 

V _ ^ i(f/c — yv) — j(x c — Xy) 

at c due to v 2a [(*„ - i,) 2 + (jfc - y,Y) 1 ' 

To extend the expression to apply to the vortices and control points distributed as 
per equations 2.16 and 2.17, simply insert the appropriate indices. Then the solution 
of the vortex lattice can be formulated into a mathematical statement. The velocity 
induced at the n th control point by all the N vortices dotted with the normal vector 
is equal to minus the free stream velocity dotted with the normal vector. 



-1 



N 



£r. 



i (</c(«) - yv{m)) — j(ar c (n) - x „(m)) 



n„ = -U in fi • n n (2.20) 



27r m=i [(*c(n) - x v (m)) 2 + ( y c {n ) - y v {m)) 2 ) 

Equation 2.20 can be written for each of the control points. This represents a set 
of N simultaneous equations which can be solved for the unknown vortex strengths 
(T n ). Separating out the T n terms as a separate component, the terms that are left 
formulate a square matrix that represents the induced velocity caused by each point 
vortex on every control point. This matrix is called the influence coefficient matrix 
(A). Equation 2.20 can be cast in matrix form as: 



a r = u 



( 2 . 21 ) 



where T is the array of point vortex strengths and U is the array of boundary condi- 
tions at the control points. The A matrix can be reduced using Gaussian elimination 
or any convenient algorithm for matrix inversion[4]. Once the matrix is inverted, the 
vortex strengths can be determined. 

r = A _1 -U (2.22) 
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This completes the treatment of the 2-D vortex lattice method for a foil in an 
unbounded fluid. The next aspect to consider is the development of wall corrections 
for when a vortex lattice is located near a wall. 

2.4.3 Method of Images Applied to VLM 

Many of the cases analyzed involve 2-D foils which are in close proximity to walls. 
This is the case for testing a foil in a water tunnel. Therefore, the unbounded vortex 
lattice method needs to be modified to account for the presence of the walls. There 
are two options for modeling the walls. The first is to put vortex lattice panels on 
the walls in a similar fashion to panelizing the foil in the unbounded method. The 
second is to use the method of images where the walls are treated as “mirrors” and 
imaginary image foils are placed outside the flow domain [28]. 

Adding panels to the walls is not an attractive option. Adequate representation of 
the walls with vortex panels requires many more panels than are already used for the 
foil. This adds significantly to the computer code solution time. Furthermore, the 
boundary condition at the walls (V • n = 0) is only satisfied at the individual control 
points, rather than continuously along the boundary. Another issue with paneling 
the walls is deciding how far up and down stream to extend the panelization scheme. 
The required distribution of wall vortex panels varies and is dependent upon many 
parameters. There is no obvious generalized scheme for implementing a rule based 
algorithm to determine the minimum required number of panels and their respective 
spacing on the walls. 

A better way to represent the walls in an inviscid flow solver is to use the method 
of images when the geometry permits. In this case, the foil is located between two 
parallel walls. Therefore the image geometry is very simple. Using the wall as a 
symmetry plane, an image body is placed on the opposite side of the wall from the 
actual body. For a vortex situated near one wall the image formulation is simple - - 
add one image vortex[28]. However with two parallel walls, the image geometry is a 
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little more complicated. The two parallel walls create an infinite number of reflective 
planes. This geometry, with the first pair of vortex lattice images, is shown in Figure 
2-5. Newman[24] provides a treatment for the similar case of a potential flow point 
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Figure 2-5: Vortex Lattice Near Walls With Images 



source situated between two walls located at y = ±|6. Implementing this for the case 
of a distribution of point vortices is analogous. This requires the addition of an infinite 
array of image vortices at y = ±6, ±26, ±36, • • •, to satisfy the boundary condition at 
the walls. A closed form solution can be formulated for the infinite array, but this 
is unnecessary since each subsequent image pair has rapidly diminishing impact on 
the solution. For all practical purposes, a converged solution can be obtained with a 
small number of image pairs. This is demonstrated in section 3.5.2. The important 
thing to remember with this solution scheme is that the image vortices are the exact 
same strength as the original vortices. Therefore, no additional unknowns are added 
to equation 2.21. The only additional burden in computation is the addition of the 
influence of the vortex images to the influence coefficient matrix (A). Lastly, using 
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the image representation eliminates the problem of deciding how far up and down 
stream to extend the walls. Using images, the wall boundary condition is exactly 
satisfied for — oo < x < oo. 



2.5 Momentum Theory for Lift and Drag Calculations 



It has been shown that it is a simple task to calculate the the lift and drag by direct 
integration on the foil surface when using numerical foil analysis tools. However, in 
water tunnel experiments this is very difficult to do. It is very time consuming and the 
level of complexity of instrumentation required to take the surface measurements is 
impractical. The normal methods for measuring lift and drag for foils in the MIT MHL 
are “bounding box” velocity measurements using laser Doppler velocimetry(LDV)[17]. 
This is an application of fluid momentum theory. Bounding box calculations are used 
in this thesis to make comparisons between numerical calculations and experimental 
measurements. Bounding box contours can also be extracted from RANS solutions 
as a check of the surface pressure integration results. 

2.5.1 Bounding Box Analysis 



Bounding box methods are a way of calculating the forces acting on 2-D hydrofoils 
expressed in terms of integrals of the velocity flow field along a closed contour sur- 
rounding the foil[20]. Figure 2-6 is an example of a closed integration contour around 
a foil. Integrating around the contour C with an outward pointing unit normal vector 
n, the momentum flux M, out of C is: 



M = 




(2.23) 



where q is the fluid flux across C and p is the fluid density. The momentum flux 
through the surface of the foil is zero, since it is by definition a rigid boundary. The 
momentum flux of the volume of the fluid V enclosed on the outside by C and on 
the inside by the foil surface is given by equation 2.23. Using Newton’s 3rd law, this 
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Figure 2-6: Rectangular Contour Around 2-D Hydrofoil 

momentum must be balanced by an opposing force F which is given by: 

F = — jf Pnds + j> rtds — Fl — Fd (2.24) 

When taking measurements around foils in the MIT MHL with the LDV apparatus, 
a rectangular contour is frequently chosen. This geometry is also convenient for 
extracting data from a RANS solver or from vortex lattice field point calculations. 
The contour is chosen such that it is sufficiently far from the foil so the flow is 
considered inviscid 2 . Therefore, all the contributions due to the shear terms (r) are 
zero and the pressure terms (P) are evaluated using the Bernoulli equation [20]. 

2.5.2 Drag by Wake Defect 

Although in principle, the direct application of equations 2.23 and 2.24 should yield 
adequate results, this is not always the case for drag measurements. A special treat- 
ment must be applied in order to get a better estimate for drag by bounding box 
methods[20]. First, some ancillary variables should be defined for the ^-component 

2 A minor exception to this is the small region where the wake passes through the downstream side of the box. 
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of the velocity in different regions of the flow along the right hand side of the contour: 



u* = u, above the wake 

u M = u + A u, in the viscous region of the wake 

u * = u, below the wake (2.25) 

The first order calculation for drag is expressed as: 

F d = pUi n f [ ( u*-u)dy (2.26) 

Jrhs of C 

This can be further simplified using equation 2.25, resulting in: 

F d = pUi n f f (Au)dy (2.27) 

Jaccross wake 

This equation is correct up to the first order of A u, as long as the contour which 

crosses the wake is sufficiently down stream such that A u Uj n / holds. In the MIT 

MHL, the optical limits of the tunnel view port window do not allow for this to happen 
when large foils(c > 30cm) are tested. Therefore, a second order correction must be 
applied to equation 2.27. Typical measurements for bounding boxes around large 
foils have wake defect velocities which range as low as 0.3Ui n / to 0.5 U{ n j. Kinnas[20] 
has shown, using equation 2.24 as a starting point, that the drag including the second 
order correction is: 

F d = pUinj [ (A u)dy -p [ (A ufdy (2.28) 

Jaccross wake Jaccross wake 

These equations for bounding boxes apply to both unbounded and bounded flows. 
Direct surface integration and bounding box analyses almost never agree exactly. In 
experiments there is the inherent uncertainty associated with taking a measurement. 
Additionally, the flow is usually complicated by some unsteady phenomena that oc- 
curs such as vortex shedding. In numerical computations there are usually small 
differences associated with the how the calculation scheme is implemented. When 
comparing experiments and computer results, the lift measurements usually agree 
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quite closely. Drag results tend to vary widely. This may well be due to the fact that 
drag, in comparison to all other flow characteristics, is a small quantity. Accordingly, 
it is a hard quantity to measure with great precision and consistency. 



2.6 Reynolds Averaged Navier Stokes Equations 

2.6.1 Derivation of the RANS Equations 



One cannot discuss the RANS equations without first reviewing the equations from 
which they originate. The Navier-Stokes Equations are the equations of motion for 
a viscous fluid. The derivation can be found in many texts, of which, Sabersky[28] 
provides a clear step by step formulation. The N-S equations stem from Newton’s 
law of motion and Newton’s law of viscous friction. The equations apply to viscous 
compressible flow with varying viscosity. The N-S equations also apply to turbulent 
flow. In vector form the equations are: 

^ = --VP + i/V 2 v+ Ji/VO + f (2.29) 

j Lyb p u 

The Mach number of the cases studied is sufficiently low such that the fluid is incom- 
pressible. Therefore the 0 term is zero. 

Considering the case of 2-D turbulent flow, the velocity of the fluid can be expressed 
as a time averaged quantity with a small fluctuating term added to it. Similarly, the 
pressure can be expressed as a time averaged pressure with a small fluctuating term. 

u = u + u' 
v = v + v' 

P = p + P' (2.30) 



The fluctuating terms u' ,v' and P' are defined such that their time averaged value 



is zero. Substituting equations 2.30 into the N-S equations, yields (written for x 



component, similar for y): 

du du! du 2 duu' du' 2 duv 

1 1 1- 2 1 1 b 

dt dt dx dx dx dy 



duv' du'v du'v' 
dy + dy + dy 
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5V 



1 dP 1 dP_ d*u 5 V 

p dx p dx ^ V dx 2 + dy 2 + ^ 5a: 2 + dy 2 



)• 



(2.31) 



If the time averaged is taken for this equation, and the continuity equation is applied 

for u and v the equation simplifies to: 

Du 1 dP 1 5 ( du — 1 5 ( du — 

~Dt ~ ~~p~fa + ~ P [ p !h ~ pu ) + ~ P ~ puv ) 



(2.32) 



This equation is exactly equation 2.29 with the addition of two terms pu' 2 and pu'v' . 
The time averaged pressure and velocity satisfy the original N-S equation. In laminar 
flow these terms are zero. These extra terms are important in turbulent flow. They 
are commonly called the turbulent stresses or Reynolds stresses. 



2.6.2 Turbulence Modeling 



With the addition of the Reynolds Stress terms to the N-S equations, more equations 
must be included in the solution to evaluate the extra unknown terms. As these terms 
represent the turbulent characteristics of the flow, the extra equations required are 
usually called turbulence models. The word model is used because the formulation 
of the equations is typically based on empirical studies of turbulent flow. There are 
several models available for use in solving the RANS equations. Three of the most 
common derived models are: 



• Prandtl mixing length “zero equation ” model 

• The “two equation” k — e model 

• Baldwin-Lomax algebraic model 

All the above methods employ the concept of eddy viscosity uj. The models use 
empirical relations to establish a value for ut which is fed into the RANS equations. 
The models used by the incompressible RANS solver in this thesis are the k — e and 
the Baldwin-Lomax Model. 

The Prandtl model uses a single parameter to determine the eddy viscosity. It 
is called the mixing length l. The mixing length is the characteristic length of the 
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mean eddy size which is much smaller than the fluids mean free path. It can also 
be thought of as the length normal to the main flow where momentum exchange is 
occurring. The Prandtl mixing model takes form: 



Uf = / 2 



'du 

Ay, 



(2.33) 



This equation directly relates the eddy viscosity to the mean flow. Therefore, no 
additional unknown terms are added to the problem formulation. The mixing length 
is usually determined by the experimental results for flows which are similar physically 
to those being modeled numerically. 

A derivation of the k — e model is presented in White[30]. This model is based on 
dissipation. This requires the addition of two equations to the equations governing 
the fluid flow. One is a turbulent kinetic energy equation and the other is the tur- 
bulent dissipation equation. Included in these equations are five empirically derived 
constants. These constants can be varied so that they can be applied to different 
types of flows such as wakes and jets. If the solution does not predict the flow rea- 
sonably, then the model and its parameters probably do not match the physics of the 
chosen flow [29]. The model is not intended for flows with high curvature. In lifting 
surface analysis, there are situations where the flow has a high curvature near the 
leading edge or un-separated flow adhering to an anti-singing trailing edge. For this 
reason, the model was not extensively used in this thesis. 

The Baldwin-Lomax algebraic model is very attractive because of its computa- 
tional efficiency and accuracy[29, 30]. The model is appropriate for bodies at modest 
angles of attack where flow separation is minimal or is not expected to occur. This 
eddy viscosity model defines an effective viscosity u e . 



= v + u T 



(2.34) 



The flow is divided into two regions. Different equations for the calculation of the 
turbulent eddy viscosity uj apply in each region. In the inner region the Prandtl 
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turbulent model(equation 2.33) is used, with the van Driest [30] damping model for 
computing the mixing length. In the outer region a different set of equations applies. 
The eddy viscosity is calculated by [29]: 

u T = KC C pF wa keF K ieb(y ) (2.35) 

In the equation, K , Cep, F wa k e , and Fxiebiy ) are a set of constants and coefficients 
calculated from the local mean flow characteristics and the distance from a boundary. 
Equations 2.32(written for both u and v), 2.33, and 2.35 formulate a set of four 
equations with four unknown quantities. They are the velocities u and v, the pressure 
P and the eddy viscosity uj. 

2.7 Boundary Layer Flows 

2.7.1 Characterization of Boundary Layers 

A boundary layer, as the name implies, is a thin layer of fluid adjacent to a boundary 
such as a foil surface or a wall. Within the boundary layer the fluid undergoes a 
transition from the free stream velocity outside the boundary layer to at rest on the 
boundary. The majority of the viscous effects are confined to boundary layers. 

To employ a RANS solver effectively, it is important to discretize a flow domain 
so that an appropriate number of fluid elements are contained within the boundary 
layer to capture the essential viscous effects. In modeling flow around lifting surfaces, 
failure to adequately capture the flow in the boundary layer typically results in an 
overestimation of lift coefficient and underestimation of drag. If the boundary layer 
is missed completely the results tend toward the inviscid potential flow results. Since 
capturing the boundary layer flow is so important to the accuracy of the RANS 
solution, it is important to be able to estimate the boundary layer characteristics 
beforehand. The initial estimates for the boundary layer dimensions are based on 
flat plate theory which are also suitable for application to lifting surfaces operating 
at moderate angles of attack. 
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Many texts provide relevant discussion of boundary layer theory and character- 
istics, of which Newman [24] and White[30] are well suited to the level of discussion 
required here. 

2.7.2 Laminar Boundary Layer Flows 



Consider the case of a flat plate of length / in a two-dimensional fluid domain with 
uniform flow U. Equation 2.29 reduces to the following: 



u 



u 



du 

dx 

dv 

dx 



du dv 
dx + dy 

du 1 dP .d 2 u d 2 u 

+ v d~y = -pi; + + V ] 

dv 1 dP r d 2 v d 2 v 

+ v a~ y = ~pI7 + "fe + a? ] 



(2.36) 



The boundary conditions that apply for this flow are that u — > U just outside the 
boundary layer and then u — > 0 at the surface through a small distance 8. Within 
this small distance S above the plate, the following characteristics are apparent: 



du /[A 
dy ~ ° ( 5 ) 
du n (U 

dx~°\ l 

.... du du 

which implies:— — 

dy dx 



(2.37) 



With u and v being 0 on the plate and in light of the relationships in 2.37, 2.36 
reduces to the following equations: 



du du 1 dP d 2 u 

U dx+ V dy p dx JrU dy 2 

n 1 dP 

0 = — x— 

p oy 



(2.38) 



The following conclusions can be drawn from equation 2.38: First, the pressure change 
across the distance 5 is not significant. Second, the pressure gradient in the x-direction 
inside the boundary layer is the same as the fluid outside[24]. In the absence of a 
pressure gradient in the x-direction, these equations can be solved numerically to 
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yield the Blasius solution for a flat plate. The concept of boundary layer thickness 
is introduced as the distance 8 where the velocity u = 0.99£/. In this formulation, 
the boundary layer thickness grows with the one-half power of the local Reynolds 
number(i?e r ). The local coordinate x is referenced from the leading edge of the plate. 

8(x) = 4.9(^y ( 2 . 39 ) 

Another quantity that can be computed is the displacement thickness or 8 M . This 
defines an effective thickening of the body due to viscous effects which corresponds 
to a lost quantity of fluid flux within the boundary layer. 

= < 2 - 40 > 

Associated with the 8 * flux reduction, there is a corresponding loss in fluid momentum 
and a momentum thickness which is calculated by: 

${x)= rw( i ~Tf) dya> °- m (wy (2 - 4i) 

To calculate the drag on the flat plate due to viscous forces, the shear stress(r xy ) at 

the plate must be determined. 

( du\ _i 

T xy (x) = P w 0-332 pU 2 Re x 2 (2.42) 

Once the shear stress is known along the plate, the total frictional drag can be com- 
puted using equation 2.13. 

All of the above results for a flat plate can be extended to a general 2-D body 
without loss of validity, provided that the radius of curvature is much larger than the 
boundary layer thickness. 

The pressure gradient in the ^-direction does affect the boundary layer. The 

boundary layer will become separated from the body at a point when the following 
conditions are met: the shear stress is zero; upstream of that point the tangential 
velocity is positive; and downstream of that point it is negative. A separated region 
is characterized by the streamlines breaking away from the body downstream. A 
separated wake is an area of relatively high vorticity and low pressure. 
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2.7.3 Turbulent Boundary Layers 



If the Reynolds number is sufficiently high, the laminar boundary layer will transition 
or change to a turbulent boundary layer at some distance aft of the leading edge of the 
foil. In the laminar boundary layer the velocity profile is smooth and regular and pre- 
dictable. At some point, the flow becomes disturbed and unstable. Transition occurs. 
After this, the boundary layer will typically remain turbulent along the remainder of 
the surface. The boundary layer remains relatively thin. The exact transition point 
can only be estimated and its location is dependent on Reynolds number and local 
roughness of the surface. Typically in experimental practice, turbulent transition is 
forced at a prescribed location near the leading edge of the foil. This is normally 
accomplished by the addition of surface irregularities such as raised bumps at the 
desired point of transition. Near the surface, inside the boundary layer, there is a 
viscous sub-layer that is small compared to the overall boundary layer thickness. The 
region between the viscous sub-layer and the outer fluid is called the turbulent core. 

The fluid flow in a turbulent boundary layer is more complicated than the laminar 
case. Mixing occurs at the interface between the free stream and the boundary 
layer. The basic momentum equations that apply to the laminar case apply to the 
turbulent boundary layer as well. At the surface the shear stress boundary condition 
still applies(r ns a Velocity gradient). The flow in the turbulent core and at the 
interface with the outer fluid is complicated. The assumption that the slope of the 
velocity profile provides an adequate value for shear stress(equation 2.42) is no longer 
valid[28]. As a result, empirical derivations of turbulent shear stress near a wall 
are typically used. One of the most common formulations is the l /7 th power law 
approximation: 



u 

u 



T-'-'O?) 



1/7 



(2.43) 



After making the substitutions u T = yjr/p, u = Ui n j at y = 5, the following equations 
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result for shear stress at the surface(r n4 ) and turbulent boundary layer thickness(J): 

1/5 



= 0.0463 



5 = 0.058 



inf 



(£) * 
1/5 



U, 



in j y 



-J 



(2.44) 

(2.45) 



where the term a is determined by the computation of the integral: 

These equations are included as they were used as a check for validating the quantities 
determined from the RANS results. A more complete discussion of turbulent layers 
and the development of the associated equations is contained in Sabersky[28]. 

2.7.4 The y + Parameter 

Finite element cell spacing near any no-slip boundary needs to be much smaller than 
the boundary layer thickness S. This is necessary to accurately capture the velocity 
profile completely through the boundary layer to the no-slip surface. The distance 
that the close-in grid points are measured from the no-slip surface are measured in 
multiples of y + where: 

y + = «r- (2.47) 

V 

and 

(2.48) 



u T = 



pU 2 



Based on Anderson[2], Black[6] and present RANS methodology in the MIT-MHL, 
y+ values for the first three cells off the no-slip surface should be < 1, < 2, and < 4 
respectively. Additionally, there should be another 10 to 15 cells spaced from y + ~ 4 
to y + 100. It should be noted that failure to adhere to this tight spacing criteria 
can result in large errors. 

2.7.5 The Law of the Wall 

The Law of the Wall is a commonly used formulation for determining the velocity pro- 
file characteristics near a wall[28]. The main assertion is that the velocity profile(u) 
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is determined by the conditions at the wall(r na ), the fluid properties^, is) and the 
distance from the wall(t/). Using the Buckingham Pi theorem, two dimensionless 
parameters can be formulated from these quantities[23]. One is a non-dimensional 
velocity and the other a non-dimensional distance. There is an implied simple func- 
tional relationship between the two parameters: 



U T \ V J 



(2.49) 



As a matter of convenience, the terms in equation 2.49 are often replaced by u + = 
f{y + ). Several approximate formulas have been deduced based on experimental data. 
Of those available, the Spalding formula is quite suitable because it fits experimental 
data well from y + « 100 all the way down to the surface[30]. 



y+ = u + + e~ KB 



— 1 — KU + — 



( KU + ) 2 (KU + ) 3 



(2.50) 



2 6 

In the above equation, values for the constants k and B are 0.41 and 5.0 respectively. 
This formula is plotted in Figure 2-7. The Spalding formula for the Law of the Wall is 




Figure 2-7: Spalding Formula for Law of the Wall 

used for validation of RANS calculations in selected cases to ensure that the velocity 
profiles close to the wall are correctly captured. One measure of the quality of a 
numerical finite element grid is how well the calculated boundary layer profiles agree 
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with the Spalding formula. It can be concluded that if the grid spacing conforms to 
the y + spacing criteria and the computed boundary layer profiles conform to the Law 
of the Wall, then the numerical grid is an adequate model for the flow. 



45 



Chapter 3 

Numerical Methods Development 



3.1 Overview 

Several 2-D foil analysis computer codes are employed in this study. The majority of 
the analysis was conducted using the RANS flow solver DTNS2D[29]. Other codes 
used are PAN2D[22] and XF0IL[11]. 

DTNS2D is a generalized 2-D domain incompressible RANS fluid flow solver devel- 
oped at the U.S. Naval Surface Warfare Center, Carderock Division. DTNS2D uses a 
finite volume finite difference formulation. In a finite difference solution method, the 
flow domain is divided into a set of discrete control volumes. Within each element, the 
governing equations of the flow domain are solved in their integral form. It is capable 
of modeling foils in bounded and unbounded domains. It should be noted that there 
are several RANS computer codes that could be applied to the cases studied here. 
However, DTNS2D was a readily available code within the MIT-MHL that required 
no modification for employment in this thesis. It has been used previously for ana- 
lyzing lifting surfaces with results that agree with experimental data[25]. Therefore, 
it was deemed adequate for this study as well. 

PAN2D and XFOIL use potential flow panel methods for calculating the inviscid 
flow characteristics around a foil in an unbounded 2-D domain. The solution is then 
coupled with an integral boundary layer method to determine the viscous properties 
of the flow[15]. The two panel methods were used to validate DTNS2D results for 
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unbounded flows. The two panel methods are currently not capable of modeling foils 
bounded by tunnel walls. 

In this chapter the development process for the FIT2D computer code that provides 
the fluid geometry input files for the DTNS2D solver is presented. This involves 
discretizing the domain into a mesh of four sided polygons finite volume elements. 
Grid generation schemes and methodology are discussed. As a subsection of FIT2D, a 
2-D vortex lattice computer code[l8] was implemented to obtain an inviscid solution 
for the lift coefficient of a foil in unbounded and bounded flows. The vortex lattice 
solution is also used to grow a dividing streamline downstream from the trailing edge 
to define a RANS zonal grid boundary. 

3.2 Overview of Grid Generation 

3.2.1 Fundamental Concepts 

The overlying governing directive of grid generation is simple: discretize the domain 
into a sufficient number of elements such that the essential physics of the flow is 
captured. The implementation of this is difficult. Figure 3-1 is a representation of a 
discretized fluid domain surrounding a hydrofoil inside a tunnel. It is important to 
notice that the grid spacing is not uniform. The size of the individual elements must 
be varied depending on the expected local characteristics of the flow. For example, 
near a wall or on the surface of the foil, the vertical spacing must be quite small 
in order to accurately capture the velocity gradient within the boundary layer. For 
hydrofoils it is also important to cluster elements near the leading edge to capture the 
strong pressure gradients associated with the flow stagnation point. At the trailing 
edge, clustering is used where separated flow is expected. Additionally, the wake 
behind the foil has steep velocity gradients that slowly dissipate with the flow. The 
clustering scheme is advantageous if the small elements are located within the wake 
region. Conversely, there are regions in the flow where very little is occurring. Velocity 
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Figure 3-1: Sample of the discretized flow domain of a foil in a water tunnel. (Only 1/3 of the total 
grid lines are reproduced in the lower frame) 



and pressure gradients are near zero. In these areas, there is no reason to have finely 
spaced cells. 

In Figure 3-1, the element dimensions vary throughout the domain. The change 
in size from largest to smallest element is approximately four orders of magnitude. 
This sample representation requires 41,238 elements. In a grid representation of an 
unbounded flow domain around a foil, the element dimensions can range as much a 
six orders of magnitude. Of course one could simply find the smallest dimensional 
requirement for an element (eg. foil leading edge) and set all elements to that size and 
capture the essential physics of the flow. The problem is, for the same sample above, 
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the grid would require 6xl0 10 elements! In a typical numerical solver, the solution 
time T is proportional to the number of elements squared ( T oc n 2 ). For the same 
results, the processor time to obtain a solution has grown by a factor of 10 12 . Clearly 
this is not a practical solution. So, there are two conflicting requirements which must 
be met in order to produce results which are both practical and accurate. The first is 
to maximize the number of elements to capture all the pertinent physics. The other 
is minimize the number of elements to obtain the minimum solution time yet also 
obtain correct results. To satisfy both of these driving forces, the element sizes must 
vary throughout the domain in order to capture the local effects and minimize the 
overall number of elements. 

3.2.2 Geometric Limitations of DTNS2D 

Before proceeding with the description of the FIT2D program, it is important to spec- 
ify the geometric limitations for the DTNS2D input files. In the DTNS2D program 
the individual elements are grouped into zones. There can be any number of zones in 
DTNS2D. Figure 3-2 is a typical zone for DTNS2D. A zone consists of four sides. The 
sides can consist of any shape to enclose the zone so long as the four sides together 
comprise a simple connected region. A simple closed region has the property that an 
arbitrary closed curve lying in the region can be shrunk continuously to a point in 
the region without passing outside of the region[14]. Each pair of opposite sides must 
have the same number of elements along the side. Adjacent sides are not required to 
have equal number of elements. The elements within each zone cannot overlap any 
other elements. Every element within the zone must have four sides of finite non-zero 
length. Each zone side has a defined boundary condition (eg. no-slip, tangent, free 
stream, continuity of velocity gradient or pressure gradient). The boundary condition 
along a side can not change. Adjacent sides can have different boundary conditions. 
An example of a DTNS2D zone model is presented in Figure 3-3. This type of model 
is frequently called an “H” grid. 



49 



Typical Wall Boundary 
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Figure 3-2: Sample geometry of a DTNS2D zone. (Only 1/3 of the total grid lines are reproduced) 

3.3 Development of the Computer Code FIT2D 

3.3.1 Functional Requirements 

Recent research projects in the MIT-MHL water tunnel have demonstrated a re- 
quirement to use a RANS computer code to model flow around hydrofoils[19]. Some 
experiments on foils with unique camber distributions yielded interesting results. It 
is anticipated that the RANS output could help provide some insight into the results 
that were being observed. 

Before the development of FIT2D, generating the DTNS2D RANS solver input files 
was time consuming and tedious. A requirement was established for the development 
of a computer program that could quickly generate the flow domain and RANS input 
files for any 2-D hydrofoil that could conceivably be tested in the MIT-MHL. This 
stated need was translated into the following set of functional requirements or program 
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Figure 3-3: Typical DTNS2D zonal boundaries for a foil in a tunnel 



objectives: 

1. Run from within a simple user interface. 

2. Support basic inputs such as angle of attack, domain extents, and locations of 
walls if any. 

3. Allow the user to specify grid spacing in key areas such as the leading edge and 
trailing edge and inside the boundary layer. 

4. Evaluate user spacing of elements inside the boundary layer and compare against 
flat plate calculations. Make a recommendation for any changes. 

5. Read in x,y foil offsets. Spline the offsets and locate the leading edge. 

6. Rotate the foil to user prescribed angle of attack. 

7. Establish DTNS zone boundaries for a six zone “H” grid. 

8. Allow for the zone boundary going downstream from the trailing edge to follow 
an established wake or a first guess wake from a vortex lattice calculation. 

9. Generate all interior grid points for the mesh using isoparametric interpolation. 

10. Write datafiles for the program TECPLOT containing a representation of the 
grid that can be viewed for correctness. 

11. Write datafiles for the INMESH grid smoothing program. 
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3.3.2 Structure of FIT2D 



FIT2D is written using the FORTRAN 77 code standard. It is executed interactively 
from a UNIX command window. The typical setup when running FIT2D is to have the 
data input file or control file displayed in an EMACS window(or any other editor), and 
a TECPLOT graphics window for viewing meshes. The three window configuration 
was chosen in favor of writing a dedicated X-Windows application. A UNIX command 
window, EMACS and TECPLOT provide an environment to which many people are 
already accustomed. After setting up the initial foil offset file(/name.foil) and control 
file(/name.ctrl), the user loops through the following basic procedure: 

1. Initialize FIT2D 

2. Follow prompts to generate a grid. 

3. View the grid in the graphics window. 

4. If the grid is satisfactory, confirm permission to write INMESH files. 

5. If the grid is unsatisfactory, make changes to fname . ctrl file and reinitialize 
FIT2D. 

Once satisfactory results have been achieved with FIT2D, the program INMESH is 
run to smooth the grid. This process is described in greater detail in section 3.6. 
Once the grid has been generated and smoothed, the geometry is ready for input to 
the RANS solver. A complete listing of the Program FIT2D is contained in Appendix 
A. 

3.3.3 Selecting Required Input 

Sample input files for FIT2D are contained in Appendix B. With the exception of 
input/output control, and boundary layer resolution changes, all FIT2D input is in 
batch form. The control file contains all of the data for describing the flow domain 
and how the mesh of elements will be distributed. The required input stems from the 
ideas presented in section 3.2.1. In the interest of making grid generation a somewhat 
more mechanical and less artistic process, it was decided early on to limit user input to 
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Table 3.1: Batch Input Variables for FIT2D 



FOILGEO 


/name. foil geometry file name, located in same directory 


AOA 


Angle of Attack in degrees 


Xpiv 


distance x/c of pivot point from leading edge 


Ypiv 


distance y/c of pivot point from nose-tail line 


TUNPARAM 


foil scale parameter c/tunnel width 


USL 


upstream limit of flow domain in chord lengths from L.E. 


DSL 


downstream limit of flow domain in chord lengths from T.E. 


PHI1 


forward rake distance of vertical zone lines from leading edge 


PHI2 


aft rake distance of vertical zone lines from trailing edge 


NUS 


number of elements in x-direction upstream 


NDS 


number of elements in x-direction downstream 


NVS 


number of elements in y-direction above and below foil 


NTOP 


number of elements along chord on top surface of foil 


NBOT 


number of elements along chord on bottom surface of foil 


RESLE 


Element width in chord lengths at leading edge 


RESMID 


Element width in chord lengths at mid chord 


RESTE 


Element width in chord lengths at trailing edge 


RESWALL 


Element height in chord lengths at the walls 


RESBL 


Element height in chord lengths on the foil surface 


PACK 


fraction of NVS elements targeted within boundary layer 


Re* 


chord based Reynolds number 


MESHFILE 


INMESH input file name 


RESTFILE 


INMESH restart file name 


NUMIT 


number of INMESH iterations 


TOL 


convergence tolerance limit for INMESH 


NRANSWK 


wake adaption flag: < 0 straight line, 0: use VLM, > 0 use RANS data 


RANSWKPTS 


x, y coordinates of the RANS wake line data 



a few key parameters to guide the program. Within the control file, there are fifteen 
essential parameters for grid generation that the user must specify. Additionally, 
there are some required administrative inputs. In the control file, each line of input is 
preceded by a descriptor line. The key inputs and their functions are summarized in 
table 3.1. Most of the parameters are intuitive. However, two parameters in particular 
need further explanation. The PHI1 and PHI2 rake angles are required in the mesh 
so that elements at the leading and trailing edge do not compress to a zero volume. 
An example of a problem leading edge is one that is circular with a large radius. If 
the zone boundary extended vertically from the leading edge, then the cells to the 
right of the boundary and above the leading edge would be skewed almost ninety 
degrees. In other words, cell volume is zero making it a singularity point. Referring 



53 



to Figures 3-1 and 3-3, one can see that the zonal boundaries are raked. The visible 
effect is that the cells retain a rectangular profile as they wrap around the leading 
edge. The upstream cells experience a minimal volume compression. But, the overall 
effect is to improve the consistency of the cell volumes. The raking capability was 
added at the trailing edge to accommodate highly curved anti-singing trailing edges. 

3.3.4 Splining and Element Distribution Methods 

All of the zonal boundary lines are splined parametrically in arc length using a cubic 
splining routine. The cubic coefficients are determined and stored for later use. Cubic 
splining in arc length was not required for any of the straight line boundaries. It was 
used anyway because it allows for future growth of the computer code to easily accept 
zonal boundary lines of arbitrary shape. An example of this is using a “C” shaped 
zone that wraps around the foil. 

The user specifies how many points are to be distributed along each boundary as 
well as the endpoint element dimension. Several schemes were explored to distribute 
the remaining points along the boundaries. The rate of change in cell dimensions 
should be relatively uniform along each boundary. Uniform spacing obviously is not 
a candidate. Cosine spacing was investigated, but it proved inadequate because the 
growth rate of cells near the areas of key interest was too fast. Ultimately, two 
spacing methods proved preferential for distributing elements. One spacing method 
is for horizontal zone lines. The other is for vertical zone lines. 

The horizontal lines define the tunnel walls, foil surface, upstream zone lines, and 
the wake. Along horizontal lines a parametric cubic distribution routine is employed. 
This routine was developed by Black[5] for distributing points on axisymmetric bodies 
and ducts. The program has been converted to a subroutine for inclusion in FIT2D. 
It has proven to be adequate for a variety of geometries. 

The spacing of elements on the vertical lines is somewhat more difficult. Vertical 
lines define how the elements are distributed normal to the foil and the tunnel walls. 
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There are two driving factors. One is to place an adequate number of very small 
elements close to the foil and walls to accurately capture the boundary layer profile. 
The other is to have much larger elements outside the boundary layer to minimize 
the total number. The user sets the spacing by specifying values for the variables: 
RESBL, RESWALL, NVS, and PACK. RESBL and RESWALL should be set using 
the y+ method of section 2.7.4. Poor selection of RESBL and RESWALL (i.e. too 
large) will yield poor results and the boundary layer characteristics will be missed. 
The next trade-off is with NVS and PACK. It is always easy to increase number of 
points NVS, but the solution time suffers. PACK dictates how many of the total 
points NVS will be allocated to the boundary layer. 

Given that the steepest velocity gradients occur nearest to the surface, the fol- 
lowing spacing routine for vertical lines has been chosen. Set first element height 
at RESBL or RESWALL. Then geometrically grow the cells by a factor of 1.25 un- 
til the number of elements corresponding to PACK are used. Lastly, distribute the 
remainder of the cells using the cubic distribution routine. 

3.3.5 Defining Domain Boundaries 

Defining the length limits of the discretized domain was approached in several ways. 
For analyzing cases for foils in the tunnel, the actual physical limits of the tunnel 
were used to constrain the upstream and downstream limits of the gridded geometry. 
In the initial formulation of FIT2D, it was not foreseen that there would be a large 
need for analyzing unbounded flow around foils using the RANS code. However, as 
the research process progressed, more and more cases were analyzed in unbounded 
domains. 

The bounded case methodology. 

Appendix C shows the configuration of the MIT-MHL water tunnel test facility. 
Figure C-2 highlights the specifics of the test section where foils are mounted for 
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measurement. Once the scale of a subject foil is known with respect to the tunnel di- 
mensions, the wall locations and upstream and downstream limits can be calculated. 
These parameters are supplied to the input files as: USL, DSL, and TUNPARAM. 
FIT2D then places all boundaries at the proper location to replicate the tunnel ge- 
ometry. This methodology proved to be adequate. Results obtained agree well with 
MIT-MHL tunnel experimental results for lift, drag and wake profiles. 

The unbounded case methodology. 

For unbounded analysis the domain limits were extended. The boundary condition 
at the walls was changed from no-slip to tangent and the walls were moved far away 
from the foil. A convergence study was conducted using the vortex lattice subroutine 
to determine how far to locate the walls from the foil such that the flow around the 
foil could be considered unbounded. A distance of five chord lengths was established 
for the walls. At this distance, the difference between the lift coefficient calculated 
by VLM with images and without is approximately 0.1%. Pressure distributions on 
a foil surface were compared against results from the computer code PAN2D and the 
differences were negligible. The flow domain was extended upstream and downstream 
as well. Five chord lengths was judged suitable. Observation of pressure and velocity 
contours from the RANS outputs showed good uniformity at the inflow and outflow 
ends of the domain. 

3.3.6 Isoparametric Interpolation 

Once the zone boundaries have been identified and the cell spacing on the boundary 
is established, the initial location for the interior cell vertices must be established. 
This is a critical step in the grid generation process. The challenge is to locate all 
the interior points such that none of the quadrilaterals overlap and all the points 
are physically interior to the zone boundary. For a simple geometric shape such as 
a rectangle or a trapezoid the task is not overly demanding. Difficulties arise when 
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adjacent boundaries are curved or highly skewed. To overcome these difficulties, a 
methodology was devised wherein the arbitrary zone boundary geometry is paramet- 
rically mapped to a unit square. The interior points are determined using simple 
line geometry. Then, the interior points are mapped back out using the inverse map- 
ping function. This process is schematically shown in Figure 3-4. Numerically, this 
is accomplished using the subroutine INTERP which is a variation of isoparametric 
interpolation adapted from Bathe’s method for finite element formulation [4]. 




3.4 FIT2D Output 

3.4.1 Tecplot Files 

Two ASCII plotting data files are generated by FIT2D. These files conform the data 
format required by the Amtec program TECPLOT version 7. 

The first file is rotate. pit. It contains the following: 
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• input geometry rotated to the desired angle of attack 

• an overlay of splined foil surface points used for foil surface zone boundaries 

• a horizontal reference line that passes through the a;, y location of the foil pivot 
point. 

The purpose of this file is to validate the input foil geometry and to ensure that the 
pivot point and angle of attack have been specified properly. A sample of a figure 
created with rotate. pit is contained in Figure 3-5. 




Figure 3-5: Display of rotate. pit data 

The second file is mesh. pit. This file contains the complete grid geometry for all of 
the zones. The purpose of this file is to allow the user to examine the grid computed 
by FIT2D to ensure that the desired spacing has been obtained and that the grid 
lines are smooth. Examples of plots generated with this file are presented in Figures 
3-1, 3-2, and 3-3. 

3.4.2 INMESH Input Files 

Once a satisfactory grid geometry has been obtained, the main output files from 
FIT2D are written. These files contain all the data required for running the computer 



58 



code INMESH(section 3.6). They are called the input and restart files. The input file 
is an ASCII data file which defines the zonal boundaries and how the cell points are 
distributed on the boundaries. The restart file is a binary formatted data file which 
contains all the interior zone point coordinates in addition to the zonal boundary 
points. 

3.5 Adjustment of Zone Boundaries in the Wake 

In section 2.4 the importance of aligning the grid zone boundaries aft of the foil to 
the convected wake was discussed. A special treatment was required to adapt the 
fluid mesh geometry to the expected wake line. Two methods were developed for 
wake adaption. The first is a vortex lattice method which finds the invicsid wake 
streamline that convects downstream from the trailing edge. A beneficial by-product 
of the vortex lattice solution is that a good approximation of the inviscid lift coefficient 
is obtained. The second method uses an initial RANS solution to extract a viscous 
wake streamline. 

3.5.1 Implementation of the ADAPT Subroutine 

A vortex lattice subroutine was implemented in the program which extracts the mean 
camber line of the subject foil in order to solve for the required circulation distribu- 
tion. Once the circulation distribution is known, the wake can be calculated. This is 
accomplished by extracting the field point velocity at the trailing edge and then incre- 
mentally marching in the direction of the field point velocity vector while extracting 
the field point velocity at each increment. This technique is visually demonstrated 
in Figure 3-6. The numerical implementation of the vortex lattice method in the 
subroutine ADAPT is a direct extension from the theory presented in section 2.4. 
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Figure 3-6: Vortex Lattice Geometry with Extracted Wake Line 



3.5.2 Adding Image Vortices 

The ADAPT subroutine was initially constructed without image vortices which would 
account for the wall effect. This proved satisfactory for analysis of unbounded flow 
around foils where there was little or no separation. However, the wake from the 
unbounded VLM solution tended to continue to convect away from the true wake 
when walls were present. Initially, it was thought that a correction could be applied 
to the unbounded VLM which would straighten out the wake. Several attempts were 
undertaken to develop schemes which redirected the velocity vector based on wall 
proximity and inviscid lift coefficient. This was made to work well for a single foil 
at a limited range of angle of attack. But the method was not robust for universal 
application and would require perturbing the parameters with each application. The 
addition of vortex images was then explored. Following the development in section 
2.4.3, pairs of image vortex lines were added to the ADAPT routine. The final 
results demonstrating the improved wake tracking are shown in Figure 3-7. The 
inset panel clearly shows how the unbounded solution continues to convect downward 
downstream while the wake of the VLM with images eventually becomes parallel to 
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the tunnel walls. 




Figure 3-7: Comparison of wakes for VLM(unbounded) to VLM with Images(bounded) 



3.5.3 Convergence of the Vortex Lattice Method 

Studying the convergence criteria of the VLM with images is a two-dimensional prob- 
lem. The goal of the convergence study is to determine a minimum number of com- 
ponents that must be modeled in the computer code in order to obtain a reasonable 
solution. First, one has to consider what is an adequate number of vortex panels to 
demonstate convergence for the lift coefficient. Next, one must consider how many 
pairs of images must be included such that the addition of more image pairs has no 
effect on the solution. The results of the convergence study are presented in Figures 
3-8 and 3-9. Based upon the trend of the two curves, the final configuration selected 
for the VLM solution employs 40 vortex panels and 16 sets of image pairs. This 
yields an error in inviscid lift coefficient of less than 0.04%. The percent error of lift 
coefficent is defined using the single precision numerically converged value for the lift 
coefficient with 640 vortex panels and 128 sets of image pairs as the reference. 

Structuring the grid zones aft of the foil using the VLM wake was useful in cases 
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Figure 3-8: Vortex Lattice Panel Convergence Curve 




Figure 3-9: Vortex Lattice Image Pair Convergence Curve 

where tunnel walls are present and there is little or no separation occurring at the 
trailing edge. The VLM wake was not a good approximation of the true wake for 
foils that had extreme camber distributions at the trailing edge or that had separated 
flow. In these cases the viscous effects on the lift coefficient become more important 
and the inviscid solution has less validity. When the viscous effects start to become 
significant, the VLM wake overshoots the true wake in the same manner that the 
VLM without images overshoots in the bounded flow case. 
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3.5.4 Adapting Wake Zone Boundaries Using RANS Output 

Despite refinements to the vortex lattice method, a significant difference persisted 
between the wake predicted by the VLM wake and the wake profile obtained using 
the RANS solution. Therefore, a method was developed to extract the wake data 
from an initial RANS solution based upon the the vortex lattice wake estimate. This 
is accomplished using the following steps: 

1. Generate grid using the VLM computed wake. 

2. Run the RANS solver for a few thousand iterations. 

3. Extract data for the wake centerline streamline using TECPLOT. 

4. Condition the data with the computer code GETWAKE. 

5. Update the /name . ctrl file with the wake data. 

6. Generate a new grid using the RANS wake data. 

7. Run the RANS solver to convergence. 

The program GETWAKE is listed in Appendix D. It is a simple conditioner that 
extracts a representative truncated set of datapoints from the TECPLOT streamline 
data. There is a major drawback to this wake correction process. It is very time 
consuming because the RANS solver must be used iteratively. It should be restated 
that this is only necessary in cases where separation has occurred near the trailing 
edge or the camber distribution on the aft part of the foil is unusual or both. A 
comparison between the two wake adaption schemes is shown in Figure 3-10. 

3.6 Using a Poisson Solver to Refine Grid 

3.6.1 The Program INMESH 

The computer code INMESH is used to refine the grid geometry output from FIT2D 
and generate the geometry input files for DTNS2D[8]. INMESH is a robust grid 
generation and refinement tool that has been used in a variety of fluid mesh geometry 
applications. INMESH uses an elliptic solver to approximate the Poisson equation 
over a prescribed domain with specified boundary values. The calculated streamlines 
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Figure 3-10: Sample Comparison of VLM Wake and RANS Wake Zone Adaption 

and equipotential lines are used to formulated the grid. The forcing term in the 
Poisson equation is used to cluster the grid near surface boundaries and to ensure 
orthogonality. 

INMESH uses the Successive-Over-Relaxation(SOR) method to solve the system 
of equations. The discretized equations are solved by iteratively sweeping through the 
domain and shifting the values at the nodes until convergence is achieved. A superb 
analogy for this process was described by Black[6]: 

The best analogy to this solution process would be to consider the initial 
guess at a grid to be a bedspring system where each node has four springs 
attaching it to its neighbors. The spacing of the nodes on the boundaries 
is controlled by the user and clustering could be seen as variable strength 
springs. The SOR method releases each node and allows it to move to its 
’natural’ position. By iteratively sweeping through the domain, the nodes 
will eventually stabilize to the final grid. 

INMESH is capable of commencing the grid generation process with only the 
boundary points specified. Then begins starts with a set of uniformly spaced interior 
points and starts the bed spring process. By employing the isoparametric spacing 
routine in FIT2D, a reasonable first guess at the location of all the interior points 
is obtained. These points are provided to INMESH in the binary restart file. A 
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significant benefit from doing this is that the INMESH code requires several thousand 
fewer iterations to achieve the same convergence tolerance. In general, a higher 
relaxation coefficient can also be used. Figure 3-11 demonstrates the results of grid 
refinement using the program INMESH. 




Ouput Mesh 
From FIT2D 
(isoparametric 
spacing) 




Figure 3-11: Typical Grid Spacing Before and After Smoothing with INMESH 



3.6.2 Correction of INMESH output for use in DTNS2D 

The INMESH code is used for a variety of grid generation applications. It provides 
the input files for the DTNS family of RANS solvers. INMESH writes out DTNS 
input files that overlap the individual zones by one cell at each adjoining boundary 
which are used if higher order boundary conditions are desired. The DTNS2D variant 
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of the code does not support this overlapping geometry configuration. Therefore a 
small computer code is used to condition the INMESH output so that it is compatible 
with the DTNS2D code. The program is called PATCH and it is listed in Appendix 
E. PATCH is specifically tailored for a six zone “H” grid scheme. 

3.7 The RANS Solver 

3.7.1 Overview of DTNS2D 

The flow solver used in this thesis is the David Taylor Navier-Stokes Two Dimensional 
computer code developed by Gorski[13], The DTNS2D computer code has been 
validated with experimental results for flow around two dimensional lifting bodies by 
Nguyen[25]. 

The DTNS2D code solves the RANS equations for incompressible turbulent flow 
which are derived in section 2.6. The flow domain is discretized into a large num- 
ber of quadrilateral cells. The solution method is a cell centered finite difference 
method based on a finite volume derivation. The program is structured so that the 
flow domain can be broken into a number of separate blocks each having its own set 
of boundary conditions. This allows for modeling of flow around foils in tunnels or 
unbounded domains. Zone junctions support a number of possible boundary condi- 
tions such as: no-slip, tangent, continuity, pressure specified, or velocity specified. A 
typical zone scheme and associated boundary conditions was demonstrated in Figure 

3-3. 

The code utilizes a third-order upwind difference total-variation-diminishing (TVD) 
scheme applied to the convection terms and a second-order central difference scheme 
applied to the diffusion terms. A first order Runge-Kutta forward time marching 
scheme is used to determine the steady state solution. 
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Chapter 4 



Data Post Processing Methods 

4.1 Introduction 

A variety of tools were used for post processing of the RANS output data. At a 
minimum, it was desired to obtain the lift and drag coefficients for a foil at a given 
Reynolds number and angle of attack. For the purposes of comparing the RANS 
solution with other foil analysis computer codes, the pressure distribution on the foil 
surface was also required. During the study, other questions arose regarding specifics 
of the flow characteristics within the boundary layer on the foil surface and in the 
wake. For studying foils in tunnels, it is beneficial to be able to extract tunnel wall 
boundary layer profiles. 

4.2 UNNS2D 

The output from DTNS2D contains all the flow characteristics at the center of each 
cell. However, the data is not in a readily viewable or workable format. A data 
conversion program UNNS2D, written by Scott Black of the MIT-MHL, was used 
to convert the data from raw output to a format compatible with TECPLOT v7.0. 
UNNS2D is a simple program and it could be easily adapted to convert data to any 
popular graphics software data format. The FORTRAN 77 source code for UNNS2D 
is located on the MIT-MHL software archive server. Once the data is loaded into 
TECPLOT, it can be displayed and manipulated in a variety of ways to view field 
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pressure and velocity contours or to draw streamlines. The UNNS2D output file is 
called dtns2d.dat. It contains data in columns by zone for each cell centered coordinate 
in the following orderrx, y, u, v and P. 

4.3 Computing Lift and Drag 

Lift and drag are computed using the program FLD2D (Foil Lift and Drag Two- 
Dimensional). It is an adaptation of a program that was previously used in the MIT- 
MHL to calculate the lift and drag for the duct and body components on ducted 
axisymmetric bodies. FLD2D extracts the cell centered flow properties from the 
DTNS2D output files. The lift and drag are numerically integrated using the methods 
of section 2.3.3. The primary outputs from FLD2D are the shear component of the 
drag coefficient and the drag and lift coefficient obtained by integration of the surface 
normal pressure. The shear component of the lift is neglected because it is a very 
small quantity when compared to the lift obtained by pressure. A plot file called 
cp.dat contains data for the pressure distribution ( Cp ) on the foil and numerical 
trail for the integrated quantities. The FORTRAN 77 source code for FLD2D is also 
located on the MIT-MHL software archive server. 

4.4 Bounding Box Calculations 

One focal point of this thesis is to evaluate differences between computational results 
and experimental measurements. In the MIT-MHL Water Tunnel the CONTOUR 
bounding box program is used to calculate the lift and drag of 2-D foils from measured 
velocities around a bounding contour. The contour integration techniques used by 
the program CONTOUR are presented in section 2.5.1. A bounding box can be 
extracted from the dtns2d.dat file using TECPLOT and subsequently processed by 
the CONTOUR program. This is useful for a variety of comparisons. For instance, a 
side by side comparison of contour normal velocities can be made between identical 
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boxes used in an experiment and an equivalent geometry RANS solution. Bad data 
points in experimental measurements can be identified. Or if good experimental 
data is available, poor modeling in the RANS code or grid scheme can be identified. 
Additionally, bounding box calculations when applied to the RANS output data, 
provide a validation of the surface pressure integration routine(FLD2D) for lift and 
drag. 

4.5 Extracting Velocity Profiles in the Boundary Layer 

In the early stages of this thesis, research focused on the validation of previously 
obtained lift and drag characteristics of foils. It became apparent that the details of 
the flow inside the boundary layer near the foil and tunnel walls were also of interest. 
Using the LDV apparatus, it is possible to take velocity data at a sufficient number of 
points inside the boundary layer to get an accurate representation of the profile. The 
same data can be extracted easily from the RANS output using the TECPLOT data 
extraction feature. Once the profiles are extracted, characteristics of the boundary 
layers can be quantified using the methods of section 2.7. 

4.5.1 Post Calculation of y + 

A validation check of the grid adequacy in the boundary layer region is obtained by 
extracting RANS velocity and position data from the first several cells above the foil 
surface. A reverse calculation of the y + value at these three locations is performed. 
The y + values for cell spacing are compared with the guidance in section 2.7.4. The 
computer code YPLUS is the small program used for this task. The FORTRAN 77 
source code for YPLUS is located on the MIT-MHL software archive server. 

The results for lift and drag are highly dependent upon the spacing of the first 
several cells above a no-slip surface. To demonstrate this a foil was analyzed using 
the FIT2D/DTNS2D analysis tool. For each case, the spacing of the first several 
cells above the foil surface was varied. The maximum y + used for the first cell was 
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approximately twenty. The minimum y + value used was less than one. All other 
grid parameters were held constant. Figure 4-1 shows plots of the pressure coefficient 
for each case. The cases corresponding to larger initial y + values yield higher lift 




Figure 4-1: Pressure Distribution Dependence on y + Parameter 



coefficients. The pressure distribution for the final case, with y + set at 0.8, compares 
well with other numerical results. 

4.6 Extracting Velocity Profiles in the Wake 

It was suspected that the turbulence models used in DTNS2D were inadequate. Data 
were extracted from the RANS solution in vertical cuts downstream of the trailing 
edge of the foil. These data were compared to identical wake cuts measured experi- 
mentally. Once again, TECPLOT was useful for extracting this data. These results 
are presented in Case Study I for the HRA foil. 
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Chapter 5 

Results 



In this chapter, two case studies are presented. The first is an analysis of a proposed 
advanced foil section at several angles of attack and a single Reynolds number. The 
flow around the foil is computed using DTNS2D in unbounded and bounded flow 
domains. Unbounded results are compared with the 2-D lifting surface analysis code 
PAN2D-BL. Bounded results are compared with experimental measurements. Differ- 
ences and similarities are highlighted. A simple correction scheme is presented which 
extrapolates bounded measurements for lift to unbounded values. 

In the second case study, a foil with a highly cambered trailing edge is analyzed at 
a single angle of attack and Reynolds number. The foil is analyzed using DTNS2D 
in unbounded and bounded flow domains. Unbounded and bounded results are com- 
pared with other computer codes and experimental measurements. This foil has 
proven to be very difficult to analyze with current numerical methods. The results 
of various computer codes are different enough that further detailed experimental 
study is warranted. The results in this case study are used to indicate regions of the 
flow around the foil where additional experimental data should be obtained. This 
case study demonstrates how RANS analysis can be useful for developing a stream- 
lined experimental test plan which will maximize the useful data obtainable while 
minimizing the costly resources and time associated with experimental study. 

These case study results are based upon DTNS2D solutions obtained from nu- 
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merical grids that demonstrated convergence with optimum y + spacing. The specific 
parameters used for the generation of these grids are presented in Appendix B. 

5.1 Validity of the Results 

5.1.1 Sources of Error 

There are several potential sources of error in the DTNS2D RANS results. There is 
no provision in the computer code for modeling the laminar regime of a boundary 
layer. DTNS2D uses a boundary layer model that starts out fully turbulent from 
the leading edge of the foil or the beginning of a wall. Along the walls this error 
is probably not significant. On the foil, this effect may or may not be significant. 
The real transition point from laminar to turbulent flow is different on the pressure 
and suction sides of the foil. The errors resulting from the omission of modeling 
the laminar boundary layer and transition are difficult to quantify. Intuitively, the 
magnitude of this error is most likely dependent upon the Reynolds number, the 
camber and thickness distributions and upon the angle of attack. 

Another phenomenon which may occur is that the turbulent boundary layer may 
revert back to a laminar flow after transition has occurred in certain regimes of 
Reynolds number. This effect is nearly impossible to quantify due to the stochastic 
nature of the process. But, nonetheless, it is another aspect of the real fluid flow that 
is omitted in the RANS solution which may have some impact on the overall results. 

5.1.2 General Comments Regarding Comparison with Tunnel Experi- 
ments 

In experimental tests, boundary layer transition on the foil is usually forced to occur 
on the foil at a specific location by the introduction of a turbulent stimulating surface 
irregularity such as bumps, rivets or a trip wire. At the location of the turbulence 
stimulators turbulent transition will occur with a fair degree of certainty. Forcing 
transition allows a more equitable comparison between experimental measurements 



72 



computer codes such as PAN2D-BL and XFOIL. By tripping the real flow in the 
experiment and specifying this same location in the computer code, one can consider 
the flows to be identical. 

Although the foil section mounted in a tunnel is two-dimensional, the flow within 
the tunnel still retains some three-dimensional effects. In a pure two-dimensional 
numerical flow solver, the side walls are completely omitted. In the real tunnel, the 
side walls tend to have uneven boundary layer growth on the areas above and below 
the foil. The magnitude of these effects are just now receiving attention and are being 
evaluated through related research in the MIT-MHL[19]. 

The RANS computer code models absolute steady state conditions that are impos- 
sible to achieve in a water tunnel. In the tunnel there are minor variations in inflow 
velocity. In the real flow unsteady vortex shedding can occur at the trailing edge. 
Very thin foils operating at high lift coefficients may twist and flex under the loading. 
The tunnel is built to standard engineering tolerances. The measured angle of attack 
in the tunnel can only be set to within several hundredths of a degree. While in the 
gridded model the geometry is precisely specified. All these issues combine and add 
to the resulting uncertainty of the experimental results. 

5.2 Case Study I: The HRA Foil 

A proposed advanced foil design created by Hydrodynamics Research Associates(HRA) 
was selected for this case study. A significant amount of experimental data has been 
recorded for this foil in various configurations in the MIT-MHL Water Tunnel test fa- 
cility. Therefore, it is a natural candidate for use in validation of the FIT2D/DTNS2D 
analysis tool. 

The HRA foil section, presented in Figure 5-1, is intended for application in pro- 
peller blade design. This foil is innovative because of its reduced camber distribution 
near the leading edge for shock free entry. The moderate camber aft allows the foil 
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to sustain a broad uniform pressure distribution on the suction side enabling it to 
operate at a substantial lift coefficient while minimizing the potential for cavitation 
inception. The foil offsets are contained in Appendix F.l. 




Figure 5-1: The HRA Foil Shape 
5.2.1 Validation Check of Grid Adequacy 

As a first check of the validity of the results obtained from DTNS2D, the y + grid 
spacing and the boundary layer profiles were analyzed. Within each set of bounded 
and unbounded runs, the only parameter varied was the angle of attack. Therefore, 
only two different FIT2D control files were necessary: one for the bounded runs and 
one for the unbounded runs. The grids were identical except for the small variations 
associated with changing the angle of attack. Table 5.1 summarizes the post calcula- 
tion of the y + parameter for for the unbounded and bounded runs. The results shown 
are for the runs conducted at AO A = —0.28°. They are representative of all angles 
of attack that were analyzed. The first cell spacing falls within the target goal and 
those remaining conform to the distribution requirements specified in Section 2.7.4. 

The next check is to ensure that the boundary layer conforms to the Law of the 
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Table 5.1: Post Processing Check of y + for HRA Foil Grid 





y + Value 


y + Value 


Cell Number 


Unbounded 


Bounded 


Target for l 5t 


1.0 


1.0 


1 (actual) 


0.6 


0.8 


2 


2.0 


2.0 


3 


3.6 


3.8 



Wall. The boundary layer profiles shown in Figure 5-2 are extracted at the mid chord 
of the suction side of the foil. These profiles are consistent over a large portion of the 
foil surface except near the leading and trailing edges. At the leading and trailing 
edges there is a significant amount of curvature and fairly steep pressure gradients. 
It is natural to expect deviation from the Spalding formula because it is based on flat 
plate theory. Based upon the y + spacing and B-L profile results, the first indication 




Figure 5-2: Comparison of DTNS2D Computed B-L Profile with Spalding Formula for HRA Foil 



is that the results should be acceptable. The next step is to compare these solutions 
with other analysis methods. 
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5.2.2 Unbounded Flow Comparison 



In the unbounded case the lift and drag coefficients obtained from the RANS calcula- 
tions are compared with PAN2D-BL. A typical comparison of the foil surface pressure 
distribution is presented in Figure 5-3. This plot is prepared for a single case of the 
HRA foil at a 1° angle of attack. Comparisons at the other angles of attack eval- 
uated are similar. Throughout a range from —2° < a < +1° the FIT2D/DTNS2D 
analysis tool agrees with PAN2D-BL to within one percent of the lift coefficient. 
FIT2D/DTNS2D slightly over-predicts the total lift which is evident from the differ- 
ence in the pressure distribution near the trailing edge as shown in the inset panel 
of Figure 5-3. Values for drag coefficient computed by RANS and PAN2D-BL are 
comparable. Drag calculated from the RANS data is about 4% higher than the 
PAN2D-BL results throughout the angles of attack studied. A summary plot of these 
results is also presented later in Figures 5-6 and 5-7. 

5.2.3 Bounded Flow Comparison 

For the bounded case the lift and drag coefficients obtained from the RANS calcu- 
lations are compared with experimental results. Throughout a range from —2° < 
a < +1° the FIT2D/DTNS2D analysis tool agrees with the experimental results to 
within a few percent of the lift coefficient. FIT2D under-predicts lift compared to 
the experimental results. Some error may be caused by position measurement error 
in the water tunnel experiment. 

Values for drag coefficient computed by RANS and from experimental measure- 
ments are not as consistent as they are for the unbounded case. Drag calculated from 
the RANS data does not follow the same trend as the experimental measurements. At 
lower lift coefficients, such as AO A < —0.5°, RANS results agree well with the exper- 
imental results. There is a large jump in the RANS computed drag for higher angles 
of attack. This is inconsistent when compared with the trends of the unbounded 
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Figure 5-3: Comparison of Pressure Distribution on HRA Foil Using DTNS2D and PAN2D-BL(Ae = 
3 x 10 6 ,AO4 = 1°) 

results and calculations. This may be caused by inadequate discretization of the flow 
domain or from a breakdown of the turbulence model under these flow conditions. A 
summary plot of these results is also presented later in Figures 5-6 and 5-7. 

5.2.4 Bounding Box Comparison 

Bounding box velocity contours extracted from the RANS data are used as input 
for the CONTOUR program. The results from the contour integration of the RANS 
data along with the surface integration results are presented in Table 5.2. The contour 
integration results compare favorably with the surface pressure integration results for 
all of the unbounded data. This indicates that the RANS solution is numerically 
consistent throughout the domain. In the bounded cases, however, the CONTOUR 
and surface pressure methods do not agree. It is not necessarily clear which method 
is in error for the bounded cases. 
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Table 5.2: CONTOUR and FLD2D Results for Lift and Drag Characteristics of the HRA Foil 



AOA 


Lift CoefF 




Drag CoefF 




Unbounded 


FLD2D 


CONTOUR 


FLD2D 


CONTOUR 


-2.00 


0.0171 


0.0173 


0.0083 


0.0098 


-0.636 


0.1623 


0.1617 


0.0083 


0.0090 


-0.280 


0.2004 


0.1995 


0.0082 


0.0090 


0.885 


0.3239 


0.3224 


0.0085 


0.0092 


Bounded 


FLD2D 


CONTOUR 


FLD2D 


CONTOUR 


-2.00 


-0.0144 


0.0168 


0.0091 


0.0090 


-0.636 


0.1628 


0.2113 


0.0089 


0.0111 


-0.280 


0.2107 


0.2630 


0.0103 


0.0103 


0.885 


0.3620 


0.3869 


0.0102 


0.0087 



A sample contour obtained from the DTNS2D results and an MIT-MHL exper- 
iment are shown in Figure 5-4. The V • n component of the fluid velocity around 
the contours are plotted for each leg of the contour. The contour normal velocities 
are nearly identical for the upstream, downstream and top contour legs. The bottom 
contour leg data does not agree aft of the mid chord of the foil. 

5.2.5 Wake Profile Comparison 

To make an observation of the wake characteristics aft of the foil, vertical cuts of 
data were taken at successive locations downstream of the trailing edge. This data 
was taken at the same locations both experimentally and from the computational 
results. These cuts are presented in Figure 5-5. The first profile is taken 0.04 chord 
lengths downstream of the trailing edge. The experimental data and the RANS data 
are nearly coincident at this location. This implies that at points very close to the 
trailing edge the RANS solver is accurately capturing the true flow conditions in the 
vicinity of the trailing edge. The successive cuts still show good agreement. However, 
it can be observed that the RANS calculated wake is diffusing at a faster rate than 
the experimental wake. This may be due to ineffective turbulence model application 
by the RANS solver in the wake region. 
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Figure 5-4: Comparison of Normal Velocity Component Around Contour Box for DTNS2D and 
HRA Experimental Results(iie = 3x 10 6 ,AOA = —0.636°) 



5.2.6 Relating Bounded Measurements to Unbounded Characteristics 

The unbounded and bounded results for lift are summarized in Figure 5-6. It is 
desirable to develop a simple scheme to relate measurements taken in the tunnel 
to their unbounded values in the absence of walls. There are several options for 
implementing such a scheme. The key decision regarding which type of scheme to use 
depends on whether you believe the experimental results or the numerical results to 
be more accurate. Rather than debate the merits and faults of the these results, it 
is assumed here, for illustrative purposes, that the experimental results are accurate 
and have been obtained with a high level of confidence. 
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Figure 5-5: Comparison of Wake Profiles for DTNS2D and HRA Experimental Results(/?e = 3 x 
10 *,AOA = -0.28°) 



Given this assumption, a method is required to correct the results obtained from 
the unbounded computer code estimates to unbounded true results. The method 
formulated here neglects side wall boundary layer effects. It is assumed the RANS 
solver is fairly accurate at capturing the magnitude of the relative effects between the 
bounded and unbounded cases. Referring back to Figure 5-6, each of the lift curves 
is linear in the regime studied. The unbounded results for DTNS2D and PAN2D-BL 
are close enough that the results can be considered equal. A least squares linear 
regression for the lift curves yields the results for lift slope and intercept presented in 
Table 5.3. Each lift curve can be represented with the following equation: 
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Figure 5-6: Summary of Results for Lift Coefficient vs. AOA for the HRA Foil 



Table 5.3: Lift Slope and Intercept Data for HRA Foil Calculations 



Method 


Slope(a) 


y-axis Intercept(6) 


MHL Expt. 


0.137 Cl! AOA 


0.267 C l @AOA = 0° 


DTNS2D Bounded 


0.131 


0.247 


PAN2D-BL/DTNS2D(unb) 


0.107 


0.234 



Cl — ox A.0 A. -1- 6, (5.1) 

where a and b are the respective slopes and intercepts. The difference in the slopes 
and intercepts of the bounded DTNS2D calculated curve and the experimentally 
measured curve represent the relative error between the estimated and actual result. 
Slope and intercept correction factors are formulated by: 



a c f = a expt — a DT N S2D (bounded.) 

b c J = bexpt — boTNS2D (bounded) ( 5 - 2 ) 
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Figure 5-7: Summary of Results for Drag Coefficient vs. AOA for the HRA Foil 



Once side wall boundary layer effects have been quantified, additional terms could be 
placed on the right hand side in equation 5.2. The correction factors are applied to 
the unbounded numerically derived lift curve in the following manner: 

( Cl,) coj . = ( a un b + a c j ) x AOA + (b un b + b c j). (5-3) 

Equations 5.2 and 5.3 are applied using the values presented in Table 5.3. Figure 5-8 
shows the corrected lift curve compared with the previously obtained results. 

What is learned from this case study is that a RANS solver can be used as a liaison 
between bounded and unbounded flows. The FIT2D/DTNS2D analysis tool agrees 
well with current state of the art unbounded foil analysis tools such as PAN2D-BL. 
In bounded flow the FIT2D/DTNS2D analysis tool agrees well with experimental 
measurements. Therefore, to conduct a complete analysis for a given lifting flow, 
where comparison with experimental results are desired, the following methodology 
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Figure 5-8: Lift Curve With Correction Factors Applied for HEA Foil 



is proposed: 

1. Based on the fact that FIT2D/DTNS2D and PAN2D-BL agree well for un- 
bounded flow, use PAN2D-BL to evaluate the foil in the range of angle of attack 
and Reynolds numbers that will be used in the experimental study. 

2. Select a few angles of attack and Reynolds numbers to evaluate using the FIT2D/DTNS2D 
tool for the bounded case. 

3. Obtain experimental results and compare the results to the bounded FIT2D /DTNS2D 
results. 

4. Using equations 5.2 and 5.3, compute the correction factors for the unbounded 
predictions and calculate the corrected lift curve. 

The above methodology should make efficient use of computational assets and result 
in a consistent comparison method between experimental and numerical results. The 
corrected lift curve provides an estimate of the error between the actual measured 
lift and drag of a foil and the numerically predicted unbounded lift and drag. The 
long term goal of the correction factor scheme is to provide motivation to improve 
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measurement and numerical prediction techniques until the correction factors are 
reduced to zero. Once the correction factors are reduced to zero then it can be 
concluded that the numerical computer codes have captured all the essential physics 
of the problem. 

5.3 Case Study II: Foil With A Cupped Trailing Edge 

The development of the “Cupped Foil” is interesting. The parent shape of the cupped 
foil is derived from a NACA 0016 thickness distribution foil with an a = 0.8 mean 
line and a maximum camber of 2.55%. The parent foil was given the name “Bl” 
and is representative of a typical destroyer propeller blade section at r/R = 0.7. 
The trailing edge was beveled using the standard U.S. Navy formula for anti-singing 
trailing edges. 

Bloch suggested that it was possible to design a foil with a substantially higher 
lift to drag ratio by shifting the maximum camber distribution aft and cupping the 
trailing edge downwards[7]. This conclusion was based on a series of numerical results 
from the computer code XFOIL. Subsequently, Jorde[16] conducted experiments to 
compare the Bl parent foil to the Bl modified with a cupped trailing edge. The 
modified Bl foil section with the cupped trailing edge is presented in Figure 5-9. The 
foil offsets are contained in Appendix F.2. 

Although the cupping at the trailing edge is not extreme to the human eye, it is 
proving to be a substantial challenge to the various array of computer codes available 
for numerical analysis. Currently there is no agreement among the computer codes 
that have been used to analyze the foil. Therefore, it is an interesting foil to use for 
this case study as a validation of the FIT2D/DTNS analysis tool. The first goal in 
this case study is to identify differences between the FIT2D/DTNS2D tool and other 
computer methods. This information will be used to provide guidance for future 
experimental research involving this foil as to where it would be productive to gather 
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detailed data. Once a suitable database of experimental data is obtained reflecting 
the actual conditions of the flow around the foil, then changes can be incorporated in 
the computer codes to make their methods more accurate. The difficulties associated 
with the unique geometry of the cupped B1 foil provide a valuable opportunity to 
improve the robustness of state of the art computational tools. 




Figure 5-9: The B-l Foil With Cupped Trailing Edge Modification 



5.3.1 Validation Check of Grid Adequacy 

As with the HRA foil, the first check of the validity of the results obtained from 
DTNS2D is to check grid spacing and boundary layer profiles. The cupped foil was 
studied at one angle of attack and Reynolds number for unbounded and bounded 
calculations. Two different control files were necessary: one for the bounded run 
and one for the unbounded run. Table 5.4 summarizes the post calculation of the 
y + parameter for for the unbounded and bounded runs. The results shown are for 
the runs conducted at AO A = 0.5° and Re = 3 x 10 6 . The first cell spacing falls 
within the target goal and those remaining conform to the distribution requirements 
specified in section 2.7.4. 
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Table 5.4: Post Processing Check of y + for B-l Cupped Foil Grid 





y + Value 


y+ Value 


Cell Number 


Unbounded 


Bounded 


Target for I st 


1.0 


1.0 


1 (actual) 


0.8 


0.6 


2 


2.1 


2.0 


3 


4.0 


3.6 



The next check is to ensure that the boundary layer conforms to the Law of the 
Wall. The boundary layer profiles shown in Figure 5-10 are extracted at the mid chord 
of the suction side of the foil. These profiles are consistent over a large portion of the 
foil surface except near the leading and trailing edges. At the leading and trailing 
edges there is a significant amount of curvature and fairly steep pressure gradients 
so it is natural to expect deviation from the Spalding formula which is based on flat 
plate theory. 




Figure 5-10: Comparison of DTNS2D Computed B-L Profile with Spalding Formula for B1 Cupped 
Foil 
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5.3.2 Unbounded Flow 



Using the post processed DTNS2D output from the unbounded model, velocity and 
pressure contours can be viewed for the flow domain. Velocity contours near the 
trailing edge are presented in Figure 5-11. On this plot flow separation can be observed 
on the suction side of the trailing edge of the foil over the last three percent of 
the chord. This separation is de-cambering the flow at the trailing edge causing a 
significant reduction in lift coefficient from the inviscid value obtained from the VLM. 
The flow is also decelerating beneath the cupped region. 





Figure 5-11: Velocity Contours Near Trailing Edge of Cupped Foil (unbounded flow) 



Pressure contours around the leading edge of the foil are presented in Figure 5-12. 
The concentric circles of pressure contour point to the leading edge stagnation point. 
An interesting observation of the pressure contours is that the value of the stagnation 
point pressure exceeds the pressure that is analytically possible. This could result 
from artificial compressibility introduced into the numerics of the RANS solver(an 
additional equation of state) to run incompressible solutions. Another possibility is 
simply numerical noise. The cell sizes near the stagnation point are very small. The 
calculations for those cells may be near the numerical precision limit of the CPU. 
This stagnation pressure error is observable in the data from all three of the RANS 
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solvers compared in this case study. 




For the analysis of the flow around the cupped foil FIT2D/DTNS2D results are 
compared with a RANS study by C.I. Yang of DTMB and results obtained by Kerwin 
using PAN2D-BL. A comparison of the calculated surface pressure coefficient for 
the three methods is presented in Figure 5-13. In general, the three methods yield 
comparable results. The largest differences occur in the leading and trailing edge 
regions. The difficulty in comparing these results relates to the difference in leading 
edge stagnation pressure. PAN2D-BL has a maximum pressure coefficient (Cp) of 
exactly 1.0, the maximum in the RANS solutions are approximately 1.1 for DTNS 
and 1.15 for C.I. Yang. Correcting for these differences may yield better alignment 
of the pressure coefficient curves. In any event, the differences in the trends at the 
leading and trailing edges are significant enough to warrant detailed experimental 
study in these regions in order to validate the predictions of the different computer 
codes. Comparison of lift and drag coefficient results are presented in Table 5.5. The 
inviscid lift coefficient obtained by the VLM is double the viscous value obtained by 
RANS and PAN2D-BL. The viscous effects on lift are significant with this foil. 
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Figure 5-13: Comparison of Pressure Distribution on Cupped Foil Obtained by Two Different RANS 
Solvers and PAN2D-BL 



Table 5.5: Comparison of Unbounded Cl & Co Calculations for B1 Cupped Foil (Re = 3 x 
10 6 ,AOA = 0.5°) 



Method 


Cl 


Cd 


VLM 


1.0555 


n/a 


FIT2D/DTNS2D 


0.579 


0.0102 


C.I. Yang RANS 


0.505 


- 


PAN2D-BL 


0.472 


- 


XFOIL 


0.5281 





5.3.3 Bounded Flow 

As is demonstrated for the unbounded flow, velocity contours near the trailing edge 
are presented in Figure 5-14. On this plot flow separation is again observed on the 
suction side of the trailing edge of the foil over the last three percent of the chord. 
In the bounded case, the separated flow region is wider and extends farther aft than 
in the unbounded case. This separation is de-cambering the flow at the trailing 
edge causing a significant reduction in lift coefficient from the value obtained from 
the VLM. As with the unbounded case, the flow is decelerating beneath the cupped 
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Figure 5-14: Velocity Contours Near Trailing Edge of Cupped Foil(bounded flow) 



Pressure contours around the leading edge of the foil are presented in Figure 5-15. 
In similar fashion to the unbounded case the concentric circles of pressure contour 
point to the leading edge stagnation point. An anomaly in the figure is the discontinu- 
ity of one of the pressure contour lines forward of the leading edge. Close inspection 
of the data file indicates the anomaly is caused by a TECPLOT interpolation er- 
ror across the zonal boundary. Additionally, as is observed in the unbounded case, 
the value of the stagnation point pressure exceeds the pressure that is analytically 
possible. 

For the analysis of the flow around the bounded cupped foil, FIT2D/DTNS2D 
results are compared with two different RANS studies by C.I. Yang of DTMB and 
Lafe Taylor of Mississippi State University. A comparison of the calculated surface 
pressure coefficient for the three methods is presented in Figure 5-16. For the bounded 
case, the FIT2D/DTNS2D results are in closer agreement with the results by Yang 
than in the unbounded case. 

Comparison of lift and drag coefficient results are presented in Table 5.6. There is 
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Figure 5-15: Pressure Contours Near Leading Edge of Cupped Foil(bounded flow) 



Table 5.6: Comparison of Bounded Cl & Co Calculations for B1 Cupped Foil (Re = 3 x 10 6 , AO A = 
0.5°) 



Method 


Cl 


Cd 


VLM 


1.2017 


n/a 


FIT2D/DTNS2D 


0.6413 


0.0117 


C.I. Yang RANS 


0.605 


- 


Lafe Taylor RANS 


0.578 


- 



a similar disparity between the inviscid results and the viscous results as is observed 
in the unbounded case. Once again, the viscous effects on lift are significant with this 
foil. 

As one would expect, the presence of the tunnel walls in close proximity to the 
foil affects the lift and drag characteristics of the foil. Figure 5-17 demonstrates the 
difference in the resulting pressure distribution on the foil. To keep the compari- 
son simple, only the FIT2D/DTNS2D results are shown. The tunnel walls cause a 
flow cambering effect to the fluid flow which increases the lift coefficient from the 
unbounded value. 
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Figure 5-16: Comparison of Pressure Distribution on Cupped Foil Obtained by Three Different 
RANS Solvers 

5.3.4 Tunnel Wall Boundary Layer Comparison 

In the bounded case, the presence of the Bl cupped foil alters the boundary layer 
growth on the upper and lower walls of the tunnel. This effect is greater on the upper 
wall. A plot of the velocity contours throughout the domain is presented in Figure 
5-18. On the upper wall, ahead of the foil, the boundary layer growth follows the 
typical trend for a flat plate. Above the foil the boundary layer thins. Aft of the foil 
the boundary layer grows rapidly. These changes in the boundary layer were initially 
observed in the RANS calculations and subsequently validated by MIT-MHL water 
tunnel LDV measurements. Boundary layer profiles for the upper tunnel wall are 
presented in Figure 5-19. The DTNS2D solution shows excellent agreement for the 
boundary layer characteristics ahead of the foil. Agreement remains good above the 
foil. In the region above and just aft of the trailing edge, the DTNS2D boundary layer 
profile deviates markedly from the experimetal results. More detailed experimental 
measurements are required in this region to identify the exact location and conditions 
along the wall where the DTNS2D solution begins to deviate from the experimental 
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Figure 5-17: Comparison of Pressure Distribution on Cupped Foil for Bounded and Unbounded 
Flow 

measurements. 

Displacement thickness(<£*) and momentum thickness ( 9 ) integrals are computed 
at several locations along the upper wall. These results are pressented in Figure 5-20. 
From the leading edge to the upstream limit of the flow domain, FIT2D/DTNS2D 
agrees with the experimental results. Downstream of the leading edge, the RANS 
results differ significanlty from the experimental measurements. This may be a result 
of poor turbulence modeling in the regions of rapid pressure changes along the wall. 

5.3.5 Separated Flow 

Separation occurs near the trailing edge of the cupped foil at the angle of attack and 
Reynolds number used in this case study. The separated flow is clearly observable in 
Figures 5-11 and 5-14. Another interesting way to observe separated flow is through 
the use of streamlines. Streamlines are obtained in the TECPLOT software package 
by following velocity vectors in the domain from a specified starting point. Figure 5- 
21 shows streamlines obtained from the DTNS2D RANS solution. Separation occurs 
over the last 3% of the foil. Figure 5-22 obtained from the C.I. Yang results is 
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Figure 5-18: u Contours for Cupped Foil in Bounded Flow 



included for comparative purposes. The separation zone in Yang’s result is similar in 
length of the zone, but it is much wider across at the trailing edge. From the global 
perspective of the total flow, this difference is minimal but it is nonetheless present. 
The difference in the two RANS results indicates that detailed experimental data is 
required to validate the flow predictions at the trailing edge. 

5.3.6 Developing An Experimental Test Plan 

The cupped foil is a challenging foil to analyze numerically. There are many in- 
consistencies in the results among the different types of numerical solutions. More 
experimental data is required to validate the different aspects of the computer codes. 
Obtaining good experimental data is a very difficult task. Additionally, precious time 
and resources associated with experimental studies preclude obtaining data “every- 
where”. Here the RANS results can be applied to help identify the key areas that 
should be measured. Two benefits are realized. First, the possibility of taking data 
in areas of the flow that are not significant is reduced. Second, areas of interest that 
may be overlooked and not measured are reduced. 

The cupped foil case study yields the following guidance for the experimentalist 
to conduct a study of the B1 Cupped Foil: 
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Figure 5-19: Upper Tunnel Wall Boundary-Layer Profiles 

• The RANS solutions for the upper wall boundary-layer growth exhibit unusual 
characteristics. Boundary layer profiles should be obtained at several locations 
along the tunnel wall to confirm this, (see section 5.3.4) 

• Separation occurs at the trailing edge of the foil. The extent of separation is not 
consistent among the RANS solvers. Detailed velocity profiles are required over 
the last 3% of the suction side of the foil with additional wake cuts down stream. 

• The velocity contours on the pressure side of the foil near the “cup” are inter- 
esting. A moderate amount of data should be obtained in this region. 



95 




Figure 5-20: Upper Tunnel Wall Boundary-Layer Displacement Thickness and Momentum Thickness 

• At the stagnation point, the pressure values obtained by the RANS solvers is non- 
physical. The stagnation point should be located experimentally and pressure 
measurements in the vicinity of the stagnation point should be obtained. 

• The different numerical computer codes do not agree on the lift and drag char- 
acteristics of the foil. Lift and drag should be measured with statistical re- 
peatability using bounding box contours. The lift and drag coefficients should 
be calculated and the V • n component on each leg of the contour should be 
compared with the identical contours in the RANS domain. 

This provides the starting point for a detailed study. Based on the current numerical 
predictions, a lot of progress can be realized with the above test plan at a single angle 
of attack and Reynolds number. Once the above tasks are complete and the computer 
codes are validated, then the study can be expanded to a range of angles of attack 
and Reynolds numbers. 
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Figure 5-21: Streamlines Near Trailing Edge of Cupped Foil(DTNS2D Solution) 




Figure 5-22: Streamlines Near Trailing Edge of the B1 Cupped Foil(C.I. Yang Solution) 
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Chapter 6 



Conclusions and Recommendations 

6.1 Conclusions 

The computer code FIT2D has been developed to rapidly generate the geometry for 
the fluid flow domain surrounding an arbitrary foil shape at a specified angle of attack 
in the MIT Marine Hydrodynamics Laboratory (MHL) water tunnel. This geometry 
is provided as input data for the RANS solver DTNS2D. A suite of software tools have 
been developed to provide post processing analysis to compare the RANS solution 
with other numerical techniques and experimental measurements. 

Through the use of the FIT2D/DTNS2D analysis tool, it has been shown that 
RANS solvers are quite useful for observing the details of the fluid foil around hydro- 
foils. The RANS results coupled with experimental data provide a good validation 
database for other numerical lifting analysis tools such as PAN2D-BL and XFOIL. 

The importance of the y + parameter should be emphasized again. The cells in the 
grid near a no-slip surface must be sized to follow the guidance in section 2.7.4 in 
order to capture the significant viscous effects within the boundary layer. Failing to 
follow this guidance will undoubtedly yield questionable RANS solutions. 

Adapting the fluid mesh zonal boundaries to the wake of the foil was somewhat 
successful. Alignment allowed better concentration of the finer discretization in the 
region of the steepest velocity changes. The wake adaption will offer significant ad- 
vantages if the turbulence modeling in the wake is improved. 
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FIT2D is robust and it can be applied to a variety of foil geometries. With minor 
modification is could be used for multicomponent systems such as yacht sail and mast 
assemblies. A version of FIT2D has been converted for use to generate grids to rep- 
resent axisymmetric ducted propulsor vechicles for use with the DTNS axissyemtric 
flow solver. 

Our knowledge of the details of two-dimensional fluid flow around hydrofoils is 
incomplete. The current lifting analysis codes do not yet achieve the results that 
minimize the performance risk in proceeding from design to production of advance 
foil section propellers without extensive model testing. It is apparant from the two 
case studies presented in Chapter 5 that improvements in propeller design computer 
codes can be achieved by first looking back to the details of two dimensional lifting 
flows. When closure is brought to the remaining details of the 2-D lifting analysis 
codes, the lessons learned can be incorporated back into the propeller design codes. 
Then advanced foil sections will be able to be incorporated in future propellers with- 
out having to undergo the monumental effort that was required in the Advanced 
Technology Demonstation for propellers at the Carderock Division, Naval Surface 
Warfare Center. 

6.2 Recommendations 

The following areas require further research in order to advance the current level of 
knowledge of two dimensional fluid flow around lifting surfaces: 

• Current turbulence models in computer codes do not work well in the wake 
region behind a foil. New models based on detailed experimental tests could do 
better at matching the rate of wake diffusion as the wake travels downstream. 
Improving the turbulence models will also improve the RANS solver’s predictions 
in the boundary layer near a no-slip surface. 
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• Foils such as the Bl cupped foil are not easily modeled by current computer 
tools. More experimental study is required for these types of foils if their ad- 
vantageous lift/drag characteristics are to be exploited in the future. The ex- 
perimental results should be used to provide validation data for improved foil 
analysis computer codes. 

• FIT2D should be modified so that it works better for generating grid geometry 
for unbounded flows. Current grid spacing algorithms cluster cells too finely at 
the outer regions of the flow domain. 

• The DTNS2D RANS solver should be modified so that boundary layers start 
out laminar and then undergo transition to a turbulent boundary layer. 

• The boundary layer growth on the side walls of a water tunnel is uneven. The 
magnitude of this effect is currently unknown. Detailed measurements are re- 
quired to quantify this effect. Additionally, it would be interesting to compare 
these measurements with a RANS solution for an equivalent three-dimensional 
domain. 

• The issue of excessive calculated stagnation pressure needs to be addressed in 
DTNS2D. It is not known if the problem results from artificial incompressibility 
within the code or if the code is not adequately imposing the specified pressure 
boundary condition at the downstream extent of the domain. 



100 



ooooooooooooo 



Appendix A 



FIT2D Progam Listing 



PROGRAM FIT2D 

c 

C Prepares all input files for program INMESH for 2D foils 

C operating in the MHL tunnel at various angles of attack. 

C 

C Reads in foil geometry file, splines offsets 

C splits upper and lower surfaces adjusts to 

C angle of attack set in control file and prepares 

C input and restart files for INMESH, 

c 

2/97 JD 

Incorporated ADAPT subroutine to use vortex lattice lifting 
line to find wake dividing line to make grid wake adaptive. 



3/97 

Updated ADAPT subroutine to include image vortices due to walls 
3/97 

Updated input and adapt to accept grid line data for wake 
from rans output . 



WRITTEN BY: JOHN DANNECKER, MIT, 1996, Last Modified: 06 FEB 97 

C**************************************** ************************ ****** 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



Program executive control file: 

i i i i i i i i i i i i i i i i m i i i i i i i i i i i i i i i i ii ii i i i i i i i i i i i i i i 



NOTE: A CARRAIGE RETURN MUST BE PLACED AT THE END 

OF ALL DATA FILES 

i i i i i i i i i i i i i i i ii ii ii ii i i i i i ii i i i i ii i ii i i i i i i i i i i i i i 



"f name . Ctrl" 

The ".Ctrl" file contains the governing information 
which specifies: 

foil file: name of .naca or .foil file 
AOA: foil angle of attack ref to free stream. 

Pivot point: point about which foil is inclined 
Scale Parameter: ratio of chord length to tunnel dim 
USL: upstream flow domain limit in chord lengths 
DSL: downstream flow domain limit in chord lengths 
PHI1: leading edge zone offset number of chord lengths 
rake zone forward at tunnel wall. 

PHI 1 : trailing edge zone offset number of chord lengths 
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c 

c 

C 

C 

C 

c 

c 

c 

c 

c 



rake zone back at tunnel wall. 

Note: PHI1, PHI2 shall always be positive and these 

parameters shall not exceed 90*/, if USL and DSL 
respectively . 

Grid Resolution Parameters. 

NUS: number horizontal points upstream 
All params NDS: number horizontal points downstream 

shall be integer NVS: number vertical points above and 

and be n*8-l below foil 

points. NTOP: number of points along top of foil. 

NBOT: number of points along bottom of foil 
RESLE: specifies smallest element size near 
the leading edge 
RESMID: midchord resolution 
RESTE: specifies smallest element size near 
the trailing edge 

RESWALL: specifies smallest vertical element 
size near a wall. 

RESBL: cell height within boundary layer 

REYNOLD: chord based reynold number 



Note: horizontal spacing along wall boundaries 

is is proportional to the spacing 
along the arc of zone boundaries 
ahead along and behind foil. 



fname.dat - desired inmesh input file name 
NUMIT: NUMBER OF INMESH ITERATIONS 

TOL: CONVERGENCE TOLERANCE FOR INMESH 

NRANSWK: -1 do staight wake line 
0 do vortex lattice 
>1 use RANS data. 



File Format: FILE HEADER (limit 40 Chars) 

fname.naca or fname.foil 
AOA PIVOTX PIVOTY 

TUNPARAM Scale Parameter (Chord/Tunnel Section) 
USL DSL PHI 1 PHI2 
NUS NDS NVS NTOP NBOT 

RESLE RESMID RESTE RESWALL RESBL, PACKFACTOR 

REYNOLD 

f name . dat 

NUMIT 

TOL 



C* ***************************************************** *********************** 
C Foil Geometry Requirements: 

C 

C The leading edge of the foil is defined as the ZERO Station. 

C After foil rotation, the Pivot point becomes the ORIGIN for 

c all further calculations. 

C 

C M fname.naca M 

C A "NACA" geometry file contains the Station, Camber 

C and Thickness parameters for a foil starting at the 

C leading edge marching aft to the trailing edge. All 

C points shall be normalized values. If no leading edge 

C radius is specified, leading edge curvature will be 

C splined based on station data. Set Leading Edge 

C Radius to 1.0 for a splined leading edge. 
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c 

C File Format: 

C 

c 

C 

C 

C 

C 

C 

C 

c 

c 

C "f name. foil" 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C File Format: 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



FILE HEADER (limit 40 Chars) 

1 (integer, this flags program to NACA type geometry) 
## Chord Length 
## Number of Stations 
Thickness(tmax) , Camber(ymax) 

Station(x/C) 1/2 Thickness (t/tmax) Camber (y/Ymax) 



Leading Edge Radius Multiplier 



A M F0IL M geometry file contains the X,Y points that 
describe the surface of the foil. Foil data points 
shall be sorted so that the file marches from the 
trailing edge lower surface to the leading edge and 
then back to the trailing edge along the upper surface. 
If the two TRAILING edge points are NOT coincident, then 
fit2d.f will connect the two points with a straight line 
and specify an additional mesh zone extending from the 
blunt trailing edge downstream. 

FILE HEADER (limit 40 Chars) 

2 (integer, this flags program to FOIL type geometry) 
## Chord Length (for normalization) 

## Number of Data Points 
X Y Training Edge Lower Surface 

X Y Leading Edge 

X Y Trailing Edge Upper Surface 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



Here is the Zoning Scheme required by INMESH and DTNS2D: 
Assemble the zones for input to transfinite interpolation. 



zone 1 I zone 2 I zone 3 

i i 

i <==========\ i 

i i 

zone 4 I zone 5 I zone 6 



Here is how an indivdual zone is specified: 



Side 1 



Side 4 
ZONE "n" 

Side 2 



Side 3 
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C*************** ****************************************************** ******* 
IMPLICIT NONE 

INTEGER NPOINTS, NFTYPE , MIDPRES , MIDSUCT ,NOSE 

INTEGER NUS , NDS , NVS, NTOP , NBOT , NTEMP , I , NI , NTEMP2 , J , NTEMP3 

INTEGER NSUCT1 ,NSUCT2 ,NPRES1 ,NPRES2,NPRESSURE,NSUCTI0N 

INTEGER NUMIT , NRESTART , NWHAT , NLINE , NZONE , M , N 

INTEGER NTOPD , NVSD , NDSD , NBOTD 

INTEGER NRANSWK , K 

REAL AOA , PIVOTX , PIVOTY, TUNPARAM, USL, DSL, PHI1, PHI2 
REAL CHORDL , FULLTHK , CAMBMX , RADLE 

REAL RESLE, RESMID, RESTE, RESWALL , RESBL , YUPWALL , YBOTWALL 
REAL TEMP2, TEMPI ,TEMP3 ,TEMP4,PACKFACT0R, REYNOLD 
REAL TOL ,E1 ,RF 

CHARACTER FIN*20, PL0TFL1*20, MESHFILE*20, RESTFILE*20 

CHARACTER F0ILGE0*20, F0ILDES*40 

CHARACTER CDUM+40, QUERY* 10, JUNKTEXT*40 

CHARACTER PR0MPT2*34 

REAL PI ,TWOPI , ZERO , ONE , HALF 

PARAMETER (PI=3 . 14159 ,TW0PI=6 . 28318 , NI=200, ZER0=0.0, 0NE=1.0) 
PARAMETER (HALF=0.5) 



C 

C ARRAYS 
C 

REAL A(3),B(3) ,C(NI) ,D(NI) ,E(NI) , XI(NI) , YI(NI) ,THCK(NI) 

REAL PRES IX (NI) .PRESIY(NI) ,PRES2X(NI) ,PRES2Y(NI) 

REAL SUCT1X (NI ) , SUCT1Y (NI ) , SUCT2X (NI ) , SUCT2Y (NI ) 

REAL PRESIXS(NI) ,PRES1YS(NI) ,PRES2XS(NI) ,PRES2YS(NI) 

REAL SUCT 1XS (NI ) , SUCT1 YS (NI ) , SUCT2XS(NI ) , SUCT2YS (NI ) 

REAL PRESSUREX(NI) , SUCTIONX(NI) 

REAL PRESSUREY (NI ) , SUCTIONY(NI) 

REAL UPLINEX(NI) ,UPLINEY(NI) ,DOWNX(NI) ,DOWNY(NI) 

REAL DOWNXS (NI ) , DOWNYS (NI ) , UPLINEXS (NI ) ,UPLINEYS (NI ) 

REAL WALL IX (NI) ,WALL1Y(NI) ,WALL2X(NI) ,WALL2Y(NI) 

REAL WALL5X(NI) ,WALL5Y(NI) ,WALL6X(NI) ,WALL6Y(NI) 

REAL WALL3X(NI) ,WALL3Y(NI) ,WALL4X(NI) ,WALL4Y(NI) 

REAL VERT IX (NI ) , VERT1Y (NI ) , VERT3X (NI ) , VERT3Y (NI ) 

REAL VERT2X(NI ) , VERT2Y (NI ) , VERT4X(NI ) , VERT4Y (NI ) 

REAL VERT5X(NI) , VERT5Y (NI) , VERT6X(NI) ,VERT6Y(NI) 

REAL VERT7X(NI) ,VERT7Y(NI) ,VERT8X(NI) ,VERT8Y(NI) 

REAL Z0NE1X(NI ,NI ) ,Z0NE1Y(NI ,NI) ,Z0NE2X(NI ,NI) ,Z0NE2Y(NI ,NI) 

REAL Z0NE3X ( NI , N I ) , Z0NE3Y (NI , NI ) , Z0NE4X (NI , NI ) , Z0NE4Y (NI , NI ) 

REAL ZONESX(NI,NI) ,Z0NE5Y(NI ,NI) ,Z0NE6X(NI ,NI) , Z0NE6Y(NI ,NI) 

REAL DZ0NE5X (NI , NI ) , DZ0NE5Y (NI , NI ) , DZ0NE6X ( NI , NI ) , DZ0NE6Y (NI , NI ) 

REAL DZ0NE3X (NI , NI ) , DZ0NE3Y (NI ,NI ) , DZ0NE2X (NI , NI ) .DZ0NE2Y (NI ,NI ) 

REAL RANSWKPTS ( 2 , NI ) 

INTEGER 10(6,8,4) 

c**************************************************************************** 

c 

c 

c 

C Open Control file f name. Ctrl: default filemanme is foil2d.ctrl 
C 

C Seek Parameters describes above, use info for later output file. 

C 

C 

C OPEN INPUT FILE 

PR0MPT2 = ’PROGRAM CONTROL FILE’//’ ( ’ // ’f oil2d . crtl ’ // ’ ) = ’ 

WRITE (* , ’ (A , $) * ) PR0MPT2 
READ (* , ’ (A) ’ ) FIN 

IF (FIN(l : 1) .EQ. ’ ’) FIN = ’ f oil2d . Ctrl ’ 

WRITE (* , ’ (A) ’ ) ’OPENING FILE: ’ // FIN 

OPEN (UNIT=4 , FILE=FIN , STATUS= ’ OLD ’ ) 

C 

C 



104 



oooooooo on on ooo 



C READ IN INPUT DATA 

READ (4, ’ (A) ’ )CDUM 
READ(4, ’ (A) ’ ) JUNKTEXT 
READ (4, ’ (A) ’ )F0ILGE0 
READ(4, ’(A) ’) JUNKTEXT 
READ(4,*)A0A, PIVOTX, PIVOTY 
READ (4, ’(A) ’) JUNKTEXT 
c 

c Tunnel param for setting wall distance 
c 

RE AD ( 4 , * ) TUNP ARAM 
READ (4, ’(A) ’) JUNKTEXT 
READ (4 , * ) USL , DSL , PHI 1 , PHI2 
READ (4, ’(A)’) JUNKTEXT 
READ(4,*)NUS, NDS , NVS, NTOP, NBOT 
READ (4, ’ (A) ’ ) JUNKTEXT 

READ(4,*)RESLE, RESMID, RESTE, RESWALL, RESBL, PACKF ACTOR 
READ (4, ’ (A) 1 ) JUNKTEXT 
READ(4,*)REYN0LD 
READ (4, ’ (A) ’ ) JUNKTEXT 
READ(4, ’ (A) ’ )MESHFILE 
READ (4, 1 (A) ’) JUNKTEXT 
READ (4 , ’ (A) ’ )RESTFILE 
READ(4, ’ (A) ’ )JUNKTEXT 
READ (4,*) NUMIT 
READ (4, ’(A)') JUNKTEXT 
READ (4,*) TOL 
READ (4, ’(A)’) JUNKTEXT 

Adding capability to read RANS wake points 
RE AD (4,*) NRANSWK 

IF (NRANSWK . EQ . ZERO ) THEN 
WRITE(* ,*) * *** ’ 

WRITE(* ,*) ’**NO RANS WAKE DATA, PROCEED WITH VLM*** ’ 
WRITE(*,*) ’ *** ’ 

ENDIF 



IF (NRANSWK.LT. ZERO) THEN 
WRITE(* , *) ’ *** > 

WRITE (* , *) ’******* *STRAIGHT WAKE SELECTED***********’ 
WRITE (* , *) ’ ****N0 VORTEX LATTICE SOLN AVAIL*********’ 
WRITE (*,*) ’ *** ’ 

ENDIF 

Test to see if there are points 
IF (NRANSWK . GT.ZERO) THEN 
WRITE (* ,*) ’ *** ’ 

WRITE(* ,*) ’*******READING IN RANS WAKE DATA*********’ 
WRITE(*,*) ’*****N0 VORTEX LATTICE SOLN AVAIL********’ 
WRITE(* , *) ’ *** ’ 

RE AD ( 4 , * ) ( (R ANS WKPTS ( 1 , k ) , R ANSWKPTS ( 2 , k ) ) , K= 1 , NRANSWK ) 
wr it e ( * , * ) ( (RANSWKPTS ( 1 , k ) , R ANSWKPTS ( 2 , k ) ) , K= 1 , NRANSWK ) 
ENDIF 

********************************************************* 
*******************END OF INPUT FILE READS ************** 
********************************************************* 

OPEN ECHO FILE TO VALIDATE INPUT FILE 

c PR0MPT2 = ’CONTROL CHECK FILE’//’ (’// ’check.dat ’//’ ) = ’ 

c WRITE (* , ’ (A ,$) ’ ) PR0HPT2 

c READ (*, ’(A) ’) FIN 
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c IF (FIN(1:1).EQ. ’ ’) FIN = ’check.dat’ 

c WRITE (* , ’ (A) ’ ) ’OPENING FILE: ’ // FIN 

c OPEN (UNIT=3,FILE=FIN,STATUS=’ UNKNOWN’) 

C 

C 

c WRITE ( 3 , * ) CDUM 

c WRITE(3 ,*)FOILGEO 

c WRITE(3,*)A0A, PIVOTX, PIVOTY 

c WRITE(3 , * )TUNPARAM 

c WRITE(3,*)USL,DSL,PHI1,PHI2 

c WRITE(3,*)NUS, NDS, NVS, NTOP, NBOT.RESLE, RESMID, RESTE, RESWALL 

C WRITE(3 , * ) RESBL , PACKF ACTOR , REYNOLD 

c CL0SE(3) 

CLOSE(4) 

C******************** ************************ *********************** 

c 

C OPEN FOIL DATA FILE FOR BASELINE GEOMETRY 
C 
c 
c 

WRITE (* , * (A) * ) ’OPENING GEOMETRY FILE: ’ // FOILGEO 

OPEN (UNIT=7 , FILE=FOILGEO , STATUS= ’ OLD ’ ) 

C READ IN INPUT DATA 

READ (7, ’ (A) ’ )FOILDES 
C 

C NFTYPE SPECIFIES FOIL TYPE 

C 2 = X , Y GEOMETRY 

C 1 = NACA TYPE FORMAT 

C 

READ (7,*) NFTYPE 
READ ( 7 , * ) CHORDL 
READ (7 , * )NPOINTS 
C 

C READ IN X,Y DATA 

C 

IF (NFTYPE. EQ. 2) THEN 
DO 5 1=1 ,NPOINTS 

READ (7,*) XI (I) , YI (I) 

5 CONTINUE 
END IF 

C 

C READ IN NACA FORMATTED DATA 

C 

IF (NFTYPE . EQ . 1 ) THEN 
READ ( 7 , * ) FULLTHK , CAMBMX 
DO 6 1=1 ,NPOINTS 

READ(7, *)XI(I) ,THCK(I) , YI(I) 

6 CONTINUE 

READ (7 , *)RADLE 
END IF 
CL0SE(7) 

C 

c 

(3************************************************************************ 

c 

C IF THE FOIL IS A NACA TYPE, IT NEEDS TO BE CONVERTED TO 

C X AND Y DATA FOR FURTHER USE IN THE PROGRAM AND FEEDING 

C IN TO INMESH. 

C 

C CHECK TO SEE IF IT IS A NACA FOIL, IF SO GO TO CONVERSION 
C SUBROUTINE. 

C 

IF (NFTYPE. EQ.l) THEN 

CALL NACACONV (NPOINTS , XI , YI ,THCK ,RADLE) 

ENDIF 
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c 

c 

C OPEN ECHO FILE TO VALIDATE FOIL GEOMETRY FILE 
C 

c PR0MPT2 = 'GEO VALIDATION FILE’//’ (’//’foilgeo.dat’//’) = ’ 

c WRITE (*,’(A,$)’) PR0MPT2 

c READ (*, ’(A) ’) FIN 

c IF (FIN(1 : 1) .EQ . ’ ’) FIN = 'foilgeo.dat' 

c WRITE (* , ’ (A) ’ ) 'OPENING FILE: ' // FIN 

c OPEN (UNIT=8,FILE=FIN,STATUS=’ UNKNOWN’) 

C 

C 

c WRITE(8,*)F0ILDES 

c WRITE(8,*)NFTYPE 

c WRITE(* ,*) 'input debugger', CHORDL.NPOINTS 

99 FORMAT (F8 . 4 , F8 . 4) 

c DO 7 1=1 ,NPOINTS 

c WRITE(8,99)XI(I) ,YI(I) 

c 7 CONTINUE 

c WRITE(8,*)CDUM 

c CLOSE(8) 

C*************** *************************************** ********** ******** 

c 

c NORMALIZE FOIL DIMENSIONS IF REQUIRED 

IF (CHORDL.NE . ONE) THEN 

CALL NORMFO IL ( CHORDL , NPO INTS , X I , YI ) 

ENDIF 

c 

c 

c************************ ************************************************ 

c 

C SEND FOIL GEOMETRY TO BE ROTATED FOR ANGLE OF ATTACK 
C 

c 

c IF ( AO A. NE. ZERO) THEN 

CALL ROTANGLE(NPOINTS , XI , YI , AOA , PI VOTX , PIVOTY) 
c ENDIF 

C 

c 

c 

c IF ( AO A. EQ. ZERO) THEN 

c WRITE(* , * ) >AOA IS ZERO, ORIGIN DEFAULTS TO LEADING EDGE. * 

c ENDIF 

C** *********** ************************************************* ********** 

c 

C PREPARE DATA FILE FOR TECPLOT TO VERIFY FOIL ROTATION 
C 

C CHECK IF USER WANTS CHECK FILE: 

WRITE(*, ’ (A,$) ’) 'CREATE A TECPLOT FILE TO VIEW ROTATION?<n>: ’ 
READ(* , ’ (A) ’ ) QUERY 
C 

C debugger 

c write(*,*) query 

c 

IF ( (QUERY . EQ . ’ Y ' ) .OR. (QUERY. EQ . 'y ' ) ) THEN 

PR0MPT2 = 'FOIL TECPLOT FILE’//’ ( ’//’rotate. pit’//’) = ’ 

WRITE (* , ’ (A , $) ’ ) PR0MPT2 
READ (* , ’ (A) ’ ) FIN 

IF (FIN(1: 1) .EQ. ’ ’) FIN = ’rotate. pit’ 

WRITE (* , ’ (A) ’ ) ’OPENING FILE: ’ // FIN 

OPEN (UNIT=9 , FILE=FIN , STATUS= ’ UNKNOWN ’ ) 

WRITE(9,*) ’VARIABLES = X,Y’ 

WRITE(9,*) ’ZONE T="Input Geometry"’ 

DO 8 1=1 ,NPOINTS 

WRITE(9,*)XI(I) ,YI(I) 
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8 CONTINUE 

WRITE (9,*) ’ZONE T= M Reference Line'” 

WRITE (9 , *)-l . 0 , zero 
WRITE(9 , * ) zero , zero 
WRITE ( 9, *) ONE, ZERO 
c CL0SE(9) 

END IF 

Q***************** ********************************************* ********** 

C SPLIT FOIL IN TO FOUR ARCS. SPLINE THE ARCS AND RETURN GRADED 

C POINT ARRAYS TO DEFINE ARCS. THIS WILL DETERMINE THE MESH DENSITIES 

C FOR EACH ZONAL AREA. 

C 

C 

C IDENTIFY SPLITTING BOUNDARIES, FOILZONE returns the integer array 
c position of the splitting points. 

C 

CALL FOILZONE(XI,NPOINTS,MIDPRES,MIDSUCT,NOSE) 

NOW TAKE THE FOIL OFFSETS AND SPLIT IN TO EIGHT ARRAYS FOR 
FEEDING IN TO THE FNSPLT SUBROUTINE. 4 each of X & Y points 

CALL ARCMAKER (XI , YI , 1 , MIDPRES , PRES IX , PRES 1Y ) 

CALL ARCMAKER ( XI , YI , MIDPRES , NOSE , PRES2X , PRES2Y ) 

CALL ARCMAKER (XI , YI ,N0SE , MIDSUCT , SUCT1X , SUCT1Y) 

CALL ARCMAKER (XI , YI , MIDSUCT , NPOINTS , SUCT2X , SUCT2Y) 

************************************************************************* 

SEND TO SPLINING ROUTINE 

Spline First Interval TE to midchord pressure side 

NTEMP is number of points want to get back from fnsplt 
NPRES1=INT( (float (NBOT)-l .0)/2 .0) 



CALL FNSPLT (MIDPRES , NPRES 1 , RESTE , RESMID , PRES IX , PRES 1Y , PRES 1XS , 
+preslYS) 



debugger 

write(*,*) ’came back from arcmaker ok, nPRESl is:’,NPRESl 
write(*,*) ’Doing Aft Pressure Side’ 

write(79 , *) ’ZONE’ 
do 4990 i=l,npresl 

write(79 , *) preslxs(i) ,preslys(i) 

4990 continue 

c end debugger 

c Now spline from midchord pressure side to LE 

c 

NPRES 2 = NPRES 1+2 

c NTEMP2 is number of points going in to fnsplt 
NTEMP2 = NOSE - MIDPRES +1 

CALL FNSPLT (NTEMP2 , NPRES2 , RESMID , RESLE , PRES2X , PRES2Y , PRES2XS , 
+PRES2YS) 

C debugger 

c write(* , *) ’Doing Forw Pressure Side’ 

c write(*,*) ’number of points GOING IN TO f nsplt : ’ ,ntemp2 

c write(*,*) ’number of points from fnsplt :’ ,npres2 

c write(79 , *) ’ZONE’ 

c do 5000 i=l,npres2 

c write(79,*) pres2xs(i) ,pres2ys(i) 
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c 5000 continue 
c end debugger 
c 

C ****************** now do other half of foil, 
c Now spline from LE to midchord suction side 

C NTEMP is number of points want to get back from fnsplt 
C 

NSUCT1=INT ( (float (NTOP ) -1 . 0 ) /2 . 0 ) 

C 

c 

NTEMP2 = MIDSUCT - N0SE+1 

CALL FNSPLT (NTEMP2 , NSUCT1 ,RESLE , RESMID , SUCT1X , SUCT1Y , SUCT1XS , 
+SUCT1YS) 

C debugger 

c 

c write(*,*) 'Doing Forw Suction Side' 

c write(*,*) 'number of points GOING IN TO f nsplt : ' ,ntemp2 

c write(*,*) 'number of points from f nsplt ,nsuctl 

c write (79,*) 'ZONE' 

c do 5010 i=l,nsuctl 

c write(79 , *) SUCTlxs(i) ,SUCTlys(i) 

c 5010 continue 
c 
c 

c Now spline from Midchord Suction Side to Trailing Edge 

c 

c 

NSUCT2 = NSUCT1+2 

c NTEMP2 is number of points going in to fnsplt 
NTEMP2 = NP0INTS - MIDSUCT+1 

CALL FNSPLT(NTEMP2 , NSUCT2 , RESMID ,RESTE , SUCT2X , SUCT2Y , SUCT2XS , 
+SUCT2YS) 

C debugger 

c write(*,*) 'Doing AFT Pressure Side' 

c write(*,*) 'number of points GOING IN TO fnsplt :' ,ntemp2 

c write(*,*) 'number of points from fnsplt :' ,nsuct2 

c write(79 , *) 'ZONE' 

c do 5020 i=l,nsuct2 

c write(79, *) SUCT2xs(i) ,SUCT2ys(i) 

c 5020 continue 
c end debugger 
c 

C******************************************************************** 

C Rejoin the pairs of arcs and make this the final foil geometry 
c 

CALL ARC JOINER (NPRES 1 , NPRES2 , NPRESSURE , PRES 1XS , PRES2XS , PRES 1YS 
+ , PRES2YS , PRESSUREX , PRESSURE Y) 

C 

CALL ARC JOINER (NSUCT 1 , NSUCT2 , NSUCTI0N , SUCT1XS , SUCT2XS , SUCT1YS 
+ , SUCT2YS , SUCTI0NX , SUCTI0NY) 
c 

c******************************************************************* 

c Append the plot file with the splined data 

c 

c 

IF ((QUERY. EQ. 'Y') .OR. (QUERY. EQ. 'y')) THEN 
WRITE(9 , *) 'ZONE T="Pressure Side”' 

DO 9 1=1, NPRESSURE 

WRITE ( 9 , * ) PRESSUREX (I) , PRESSURE Y ( I ) 

9 CONTINUE 

write(9,*) 'ZONE T= M Suction Side"' 

DO 10 1=1 ,NSUCTI0N 

WRITE(9 , *) SUCTIONX(I) ,SUCTIONY(I) 

10 CONTINUE 
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CLDSE(9) 

ENDIF 

c Define upstream geometry of tunnel: 
c 

c First do horizontal line extending from leading 

c edge to forward end of tunnel. 

NTEMP3=10 

c write(*,*) ntemp3 

CALL TUNB(USL,SUCTIONX( 1) ,SUCTIONY(l) ,UPLINEX,UPLINEY,NTEMP3) 

CALL FNSPLT(NTEMP3 ,NUS , RESLE , RESMID ,UPLINEX ,UPLINEY , UPLINEXS , 
+UPLINEYS) 

Q#* ********************************************************* ************* 

c Define downstream geometry of tunnel: 

c First do horizontal line extending downstream from 

c trailing edge to end of tunnel 

c 

c 

c Adapt the GRID boundaries to the wake based on VORTEX LATTICE 
c SOLUTION, 

c 

IF(NRANSWK. EQ . ZERO) THEN 

CALL ADAPT (PRESSUREX , PRESSUREY , SUCTIONX ,SUCTIONY ,D0WNX , DOWNY , 

+ NBOT, NTOP , DSL , NTEMP3 , TUNPAR AM ) 

ENDIF 

c 

cc 

c 

c 

c if ranswk greater than 0 then there is rans data 

c 

IF THE DATA IS BAD THEN NRANSWK WILL BE SET TO -1 AND 
A STRAIGHT WAKE WILL BE USED. 

IF (NRANSWK. GT. ZERO) THEN 

CALL RANSWAKE (DSL , PRESSUREX ( 1 ) , PRESSUREY ( 1 ) ,RANSWKPTS, 

+ D OWNX , DOWNY , NTEMP3 , NRANSWK ) 

ENDIF 



If NRANSWK is set to -1 then it will do a straight wake, 
c 
c 

IF (NRANSWK.LT. ZERO) THEN 
write (*,*) >***> 

write(*,*) ’♦♦♦♦♦DEFAULTING TO STRAIGHT MAKE********’ 

WRITE(* ♦) ’ 1 

CALL TUNB (DSL , PRESSUREX ( 1 ) , PRESSUREY ( 1 ) , DOWNX , DOWNY , NTEMP3 ) 
ENDIF 
c 

c (NIN,N0UT,DS1,DS2,XI,YI,X0,Y0) 

c Spline downstream line. 

CALL FNSPLT (NTEMP3 ,NDS , RESTE , RESMID , DOWNX , DOWNY , DOWNXS ,DOWNYS ) 
c 
c 

c debugger 

c 

c write(79 , ♦) ’ZONE’ 

c do 5030 i=l,nus 

c write(79 , ♦) uplinexs(i) ,uplineys(i) 

c 5030 continue 
c 



no 



c write(79 , *) 'ZONE' 

c do 5040 i=l,nds 

c write(79, *) downxs(i) , downy s(i) 

c 5040 continue 
c 

c********************************* ************************** ************* 
c Locate the upper and lower walls of the tunnel, 
c 

CALL WALLFIND (TUNPARAM , YUPWALL , YBOTWALL) 

C 

c***** ********************************************************** ******** 

c Assign XY Points to all sections of the walls 
c 

C XBEG , XEND , YOTALL , NPO INTS , XIN , YIN , XOUT , YOUT 

c»»»»»> 
c forward upper 
c 

TEMP2=SUCTI0NX ( 1 ) -USL 
TEMP 1=SUCTI0NX ( 1 ) -PHI 1 

CALL TUNC (TEMP 1 , TEMP2 , YUPWALL , NUS , UPLINEXS , UPLINEYS , WALL1X , WALL1Y) 
c 

c forward lower 

c 

CALL TUNC (TEMP 1 , TEMP2 , YBOTWALL , NUS , UPLINEXS , UPLINEYS , 

+WALL2X , WALL2Y ) 
c»>»»»»>»> 
c aft upper 

c 

TEMP 1=SUCTI0NX (NTOP ) +PHI2 
TEMP2=SUCTI0NX(NT0P)+DSL 

CALL TUNC ( TEMP 1 , TEMP 2 , YUPWALL , NDS , DOWNXS , DOWNYS , WALL5X , WALL5Y ) 
c 

c aft lower 

c 

CALL TUNC (TEMP 1 , TEMP2 , YBOTWALL , NDS , DOWNXS , DOWNYS , 

+WALL6X , WALL6Y ) 
c»»>»»»»» 
c above foil 

c 

TEMP 1=SUCTI0NX ( 1 ) -PHI 1 
TEMP2=PRESSUREX(1)+PHI2 

CALL TUNC (TEMP 1 , TEMP 2 , YUPWALL , NTOP , SUCTIONX , SUCTIONY , 

+WALL3X , WALL3Y ) 

c 

c below foil 

c 

CALL TUNC (TEMP2 , TEMP 1 , YBOTWALL , NBOT , PRESSUREX , PRESSUREY , 

+WALL4X , WALL4Y ) 
c 
c 

c** *********************************************** ********************* 

c 

c It is very important to have the proper vertical spacing of 
c cells within the boundary layer. Determine the y+ and compare 

c with the user RESBL and RESWALL. 

CALL FINDYPLUS (RESBL , RESWALL , REYNOLD , TUNPARAM , USL , DSL) 
c 

Q* ************************************************************ ********* 

c Develop upstream vertical point spacing. Each line will consist 
c of three segments. Very fine spacing at the ends as per 
c RESWALL RESBL 

c 

C ******+D0 all the above lines first********* 

C 

c Upstream vertical line (above) 
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TEMP1=UPLINEXS(NUS) 

TEHP2=UPLINEYS (NUS ) 

TEMP3=WALL1X(NUS) 

TEMP4=WALL1Y (NUS) 

CALL TUND (TEMPI, TEMP2 , TEMP3 , TEMP4 , RESBL , RESWALL , 

+ VERT IX , VERT1Y ,NVS , PACKFACTOR) 
c 

c Leading Edge vertical (above) 
c 

TEMP 1=UPLINEXS ( 1 ) 

TEMP2=UPLINEYS( 1) 

TEMP3=WALL1X(1) 

TEMP4=WALL1 Y ( 1 ) 

CALL TUND (TEMP 1 , TEMP2 , TEMP3 , TEMP4 , RESBL , RESWALL , 
+VERT3X , VERT3Y , NVS , PACKFACTOR) 
c 

c Trailing Edge vertical (above) 

c 

TEMPl=downXS(l) 

TEMP2=downYS ( 1 ) 

TEMP3=WALL5X ( 1 ) 

TEMP4=WALL5Y ( 1 ) 

CALL TUND (TEMPI , TEMP2 , TEMP3 , TEMP4 , RESBL , RESWALL , 
+VERT5X , VERT5Y , NVS , PACKFACTOR ) 
c 

c Downstream Vertical Line (above) 

TEMP l=downXS (nds ) 

TEMP2=downYS (nds ) 

TEMP3=WALL5X (nds ) 

TEMP4=WALL5Y (nds) 

CALL TUND (TEMP 1 , TEMP2 , TEMP3 , TEMP4 , RESBL , RESWALL , 
+VERT7X , VERT7Y , NVS , PACKFACTOR) 

***********************DO ALL THE BELOW LINES************* 



Upstream Vertical Line (below) 

TEMP 1=UPLINEXS (NUS ) 

TEMP2=UPLINEYS (NUS ) 

TEMP3=WALL2X(NUS) 

TEMP4=WALL2Y (NUS) 

CALL TUND (TEMPI, TEMP2 , TEMP3 , TEMP4 , RESBL , RESWALL , 
+VERT2X , VERT2Y , NVS , PACKFACTOR) 



c Leading Edge vertical (below) 
c 

TEMP1=UPLINEXS( 1) 

TEMP2=UPLINEYS ( 1 ) 

TEMP3=WALL2X ( i ) 

TEMP4=WALL2Y ( 1 ) 

CALL TUND (TEMP 1 , TEMP2 , TEMP3 , TEMP 4 , RESBL , RESWALL , 
+VERT4X, VERT4Y , NVS , PACKFACTOR) 



c Trailing Edge vertical (below) 
c 

TEMP l=downXS ( 1 ) 

TEMP2=downYS ( 1 ) 

TEMP3=WALL6X( 1 ) 

TEMP4=WALL6Y ( 1 ) 

CALL TUND (TEMPI , TEMP2 , TEMP 3 , TEMP4 , RESBL , RESWALL , 
+VERT6X , VERT6Y , NVS , PACKFACTOR) 
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Downstream Vertical Line (bleow) 



c 
c 
c 
c 

TEMP l=downXS (nds ) 

TEMP2=downYS(nds) 

TEMP3=WALL6X ( nds ) 

TEMP4=WALL6Y(nds) 

CALL TUND ( TEMP 1 , TEMP2 , TEMP3 , TEMP4 , RESBL , RESWALL , 

+VERT8X , VERT8Y , NVS , PACKFACTOR) 
c 
c 

c ***4^****************************************************************** 



C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

cc 

c 

c 

c 

c 

c 

c 

c 



Assemble the zones AND DO isoparametric interpolation. 
Here is the zone geometry: 



zone 1 | zone 2 I zone 3 

I I 

I <=— =====\ | 

I I 

zone 4 | zone 5 I zone 6 



Here is how an indivdual zone is specified: 



Side 1 



Side 4 
ZONE "n" 

Side 2 



Side 3 



/\ 

| Array reference to enter finite 

| interpolation scheme. 

J I 

I > 

I 



ZONE ONE 

CALL BORDERS (ZONE IX , Z0NE1Y , VERT1X , VERT1Y, UPLINEXS ,UPLINEYS , 
+ VERT3X , VERT3Y , WALL1X , WALL1Y , NUS , NVS , NI ) 

ZONE TWO 



debugger 
do 5050 i=l,nvs 

write(81 ,*) vert5x(i) , vert5y(i) 
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c 5050 continue 



c end debugger 

CALL BORDERS ( Z0NE2X , Z0NE2Y , VERT 3 X , VERT3Y , SUCTIONX , SUCTIONY 
+ , VERT5X , VERT 5 Y , WALL3X , WALL3Y , NTOP , NVS , NI ) 



ZONE THREE 



CALL BORDERS (Z0NE3X, Z0NE3Y , VERT5X , VERT5Y , DOWNXS , DOWNYS , 
+ VERT7X , VERT7Y , WALL5X , WALL5Y, NDS , NVS , NI ) 



ZONE FOUR 

CALL BORDERS (Z0NE4X , Z0NE4Y , VERT2X , VERT2Y , WALL2X , WALL2Y , 
+ VERT4X , VERT4Y , UPLINEXS , UPLINEYS , NUS , NVS , NI ) 



ZONE FIVE 



CALL BORDERS ( Z0NE5X , Z0NE5Y , VERT4X , VERT4Y , wall4X , wall4Y 
+ , VERT6X , VERT6Y , pressureX , pressureY , NBOT ,NVS , NI ) 



ZONE SIX 

CALL BORDERS (Z0NE6X , Z0NE6Y , VERT6X , VERT6Y , WALL6X , WALL6Y , 

+ VERT8X , VERT8Y , DOWNXS , DOWNYS , NDS , NVS , NI ) 

C* ************************ ***************************************** 

C Open Tecplot File for viewing Mesh, 

c 

PR0MPT2 = 'Mesh View TECPLOT FILE ' // ' ('// 'mesh. pit '// ' ) = * 
WRITE (* , ' (A , $) ' ) PR0MPT2,':' 

READ (*, '(A) ') FIN 

IF (FIN( 1 : 1) . EQ . ' ') FIN = 'mesh. pit' 

WRITE (* , ' (A) ' ) 'OPENING FILE: ' // FIN 

OPEN (UNIT=82 , FILE=FIN , STATUS= 'UNKNOWN ' ) 

c topzones 

write(82,*) 'VARIABLES =X,Y' 

write(82,*) 'ZONE T= M Z0NE 1” I=',Nus,' J=',NVS,' F=P0INT' 
do 6000 j=l,nvs 
do 6010 i=l,nus 

write (82,*) zonelx(i, j) ,zonely(i, j) 

6010 continue 
6000 continue 

write(82,*) 'ZONE T= M Z0NE 2” I=',Ntop,' J=',NVS,' F=P0INT' 
do 6020 j=l,nvs 
do 6030 i=l,ntop 

write(82 ,*) zone2x(i , j ) ,zone2y (i, j ) 

6030 continue 

6020 continue 

write(82, *) 'ZONE T= M Z0NE 3’* I=',Nds,' J=',NVS,' F=P0INT' 
do 6040 j=l,nvs 
do 6050 i=l ,nds 

write(82, * ) zone3x(i, j) ,zone3y(i, j) 

6050 continue 

6040 continue 
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c 

c 



bottom xones 



write(82,*) 'ZONE T= M ZONE 4" I=’,Nus,’ J=’,NVS,’ F=P0INT’ 
do 6060 j=l,nvs 
do 6070 i=l,nus 

write (82, *) zone4x(i, j ) ,zone4y(i, j ) 

6070 continue 
6060 continue 

write (82, *)’ ZONE T="Z0NE 5" I=’,Nbot,’ J=’,NVS,’ F=P0INT’ 
do 6080 j=l,nvs 
do 6090 i=l,nbot 

write (82,*) zone5x(i, j) ,zone5y(i, j) 

6090 continue 

6080 continue 

write(82,*) ’ZONE T="Z0NE 6 M I=’,Nds,’ J=’,NVS,’ F=P0INT’ 
do 7000 j=l,nvs 
do 7010 i=l,nds 

writ e ( 82 , * ) zone6x ( i , j ) , zone6y ( i , j ) 

7010 continue 

7000 continue 
close(82) 
c 

c Generate INMESH input file: 
c 

cC CHECK IF USER WANTS CHECK FILE: 

c 

c 

wr it e ( * , * ) ’ ############################################## ’ 
write(*,*) ’ 9 

write(*,*) ’View the grid on TECPL0T, Check boundaries for’ 
write (*,*) ’wrapping. If the grid is not correct or if’ 
write (*,*) ’you want to make changes, you may abort gen-’ 
write(*,*) ’eration of the INMESH file.’ 
write(* , *) ’ ’ 

wr it e ( * , * ) ’ ############################################## ’ 
write(* , *) ’ ’ 

WRITE(* , ’ (A , $) ’ ) ’DO YOU WANT TO WRITE THE INMESH FILE?<n>:’ 
READ(* , ’ (A) ’ ) QUERY 

debugger 

write(*,*) query 

IF ((QUERY. EQ. ’Y’) .OR. (QUERY. EQ. ’y’)) THEN 



WRITE (* , ’ (A) ’ ) ’WRITING FILE: ’ // MESHFILE 

OPEN (UNIT=80 , FILE=MESHFILE , STATUS^ ’ UNKNOWN ’ ) 



FORMAT FIRST LINE 

El is the orthogonality at boundry 
El=l . 0 

Relaxation parameter 
RF=0 . 5 

Tell INMESH it is a restart 
NRESTART=2 
c 

c Unknown parameter 
NWHAT=1 
C 
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C WRITE FIRST LINE OF THE INMESH FILE 
98 format (f 4. 2, f 12.8, i7,f 8.4, 2i4) 

97 format (i3,7i3) 

96 f ormat ( i3 , 2i5 , a20 ) 
if (numit . eq. l)then 
nrestart=l 
endif 

write(80,98) El ,T0L, NUMIT, RF,NRESTART,NWHAT 
c 
c 

c * ********** ***********DQ THE FIRST THREE ZONES************** 
C 

c ***********y r i te out Boundaries for Zone 1 
c 

write(80,96) l,nus ,nvs , 9 ***Zone l*** 1 
c 

c Zone 1 Line 1 

NLINE=1 

DATA (IC(1,I,1),I=1,8)/1,5,0 > 0,1,1,0,0/ 
write(80 , 97) (IC( 1 , I ,NLINE) ,1=1,8) 
c 

call zonelines (zonelx,zonely ,nus ,nvs ,NLINE,80) 

C 

c Zone 1 Line 2 

NLINE=2 

DATA (IC( 1,1, 2), 1=1 ,8)72,1,0, 0,1, 2, 0,0/ 
write(80 ,97) (IC(l ,1 ,NLINE) ,1=1,8) 
c 

call zonelines (zonelx,zonely ,nus ,nvs,NLINE,80) 
c 

c Zone 1 Line 3 

NLINE=3 

DATA (IC(l,I,3),I=l,8)/3,0,2,l,l,l,0,0/ 
write(80 , 97) (IC( 1 , I ,NLINE) ,1=1,8) 
c 

c Zone 1 Line 4 

NLINE=4 

DATA (IC(1,I,4),I=1,8)/4,6,0,0,1,4,0,0/ 
write(80,97) (IC(1 , I ,NLINE) ,1=1 ,8) 
c 

call zonelines (zonelx,zonely ,nus ,nvs ,NLINE,80) 

Cc 

c 

c ********2one 2 

c Write out Boundaries for Zone 2 
c 

write(80,96)2,ntop+l ,nvs , 9 ***Zone 2***’ 

c Zone 2 Line 1 
NLINE=1 

DATA (IC(2 , I , 1) , 1=1 ,8)/l,0, 1,3,1, 1,0,0/ 
write(80 , 97 ) (IC(2 , I ,NLINE) ,1=1,8) 



c Zone 2 Line 2 

NLINE=2 

DATA (IC(2,I,2),I=l,8)/2,6,0,0,2,2,0,0/ 
write(80 , 97 ) (IC(2, I,NLINE) ,1=1,8) 
c 

write (80 , *)zonelx(nus-l , 1) ,zonely (nus-1 , 1) 
call zonelines (zone2x,zone2y ,ntop,nvs ,NLINE,80) 
c 

c Zone 2 Line 3 

NLINE=3 

DATA (IC(2, I, 3), 1=1, 8)/3,0,3,l,l, 1,0,0/ 
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c 

c 



write(80 , 97 ) (IC(2, I ,NLINE) ,1=1,8) 



Zone 2 Line 4 
NLINE=4 

DATA (IC(2,I,4) ,I=l,8)/4,6,0,0,2,4,0,0/ 
write(80 ,97 ) (IC(2,I ,NLINE) ,1=1,8) 

write (80, *)zonelx(nus-l ,nvs) ,zonely(nus-l ,nvs) 
call zonelines (zone2x,zone2y ,ntop,nvs ,NLINE, 80) 

c 

c 

c c** ******Zone 3 

c Write out Boundaries for Zone 3 
c 

write(80,96)3,nds+l ,nvs , ’ ***Zone 3***’ 



c Zone 3 Line 1 

NLINE=1 

DATA (IC(3,I,l),I=l,8)/l,0,2,3,l,l,0,0/ 
write(80 , 97 ) (IC(3,I,NLINE) ,1=1,8) 
c 
C 

c Zone 3 Line 2 

NLINE=2 

DATA (IC(3,I,2),I=1,8)/2,1,0,0,3,2,0,0/ 
write(80 , 97 ) (IC(3 , I , NLINE) , 1=1 , 8) 
c 

wr it e (80, *)zone2x (ntop-1 , 1) ,zone2y (ntop-1 , 1) 
call zonelines (zone3x,zone3y , nds ,nvs , NLINE, 80) 
c 

c Zone 3 Line 3 

NLINE=3 

DATA (IC(3,I,3),I=l,8)/3,2,0,0,3,3,0,0/ 
write(80,97) (IC(3, I , NLINE) ,1=1,8) 
c 

call zonelines (zone3x,zone3y, nds ,nvs, NLINE, 80) 
c 

c Zone 3 Line 4 

NLINE=4 

DATA (IC(3,I,4),I=l,8)/4,5,0,0,3,4,0,0/ 
write(80,97) (IC(3, I , NLINE) ,1=1,8) 



write (80, *)zone2x (ntop-1, nvs) ,zone2y(nus-l ,nvs) 
call zonelines (zone3x,zone3y, nds , nvs, NLINE, 80) 

ccmmmnmmm 

c 

(>************dq the second three zones******************************* 



c 



c ***********Write out Boundaries for Zone 4 
c 

write(80, 96)4, nus , nvs , ’ ***Zone 4***’ 
c 

c Zone 4 Line 1 

NLINE=1 

DATA (IC(4,I,1),I=1,8)/1,5,0,0,4,1,0,0/ 
write(80 , 97 ) (IC(4, I , NLINE) ,1=1,8) 
c 

call zonelines (zone4x,zone4y,nus , nvs , NLINE, 80) 
C 

c Zone 4 Line 2 

NLINE=2 

DATA (IC(4,I,2) ,I=l,8)/2,6,0,0,4,2,0,0/ 
write(80,97) (IC(4 , I , NLINE) ,1=1,8) 
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c 

call zonelines (zone4x,zone4y ,nus ,nvs,NLINE,80) 
c 

c Zone 4 Line 3 

NLINE=3 

DATA (IC(4,I,3),I=1,8)/3,0,5,1,4,1,0,0/ 
write(80 ,97) (IC(4, I ,NLINE) ,1=1,8) 
c 

c Zone 4 Line 4 

NLINE=4 

DATA (IC(4,I,4) , 1=1,8) /4, 1,0, 0,4, 4, 0 , 0 / 
write(80,97) (IC(4, I .NLINE) ,1=1,8) 
c 

call zonelines (zone4x , zone4y , nus , nvs , NLINE , 80) 
Cc 
c 

c ********£one 5 

c Write out Boundaries for Zone 5 
c 

write(80 ,96)5, nbot+1 , nvs , * ***Zone 5*** * 

c Zone 5 Line 1 

NLINE=1 

DATA (IC(5,I,1) , 1=1, 8 ) /l, 0,4, 3, 4, 1,0,0/ 
write(80 , 97) (IC(5 , I , NLINE) ,1=1,8) 
c 
C 

c Zone 5 Line 2 

NLINE=2 

DATA (IC(5,I,2) ,I=l,8)/2,6,0,0,5,2,0,0/ 

write(80 ,97) (IC(5, I, NLINE) ,1=1,8) 
c 

write(80,*) zone4x(nus-l , 1) ,zone4y (nus-1 , 1 ) 
call zonelines (zone5x,zone5y ,ntop,nvs , NLINE, 80) 
c 

c Zone 5 Line 3 

NLINE=3 

DATA (IC(5, I, 3), 1=1, 8 ) /3, 0,6, 1,4, 1,0,0/ 
write(80 ,97) (IC(5, I, NLINE) ,1=1,8) 
c 

c Zone 5 Line 4 

NLINE=4 

DATA (IC(5,I,4),I=l,8)/4,6,0,0,5,4,0,0/ 
write(80 , 97 ) (IC(5, I , NLINE) ,1=1,8) 
c 

write (80 ,*) zone4x(nus-l ,nvs) ,zone4y(nus-l ,nvs) 
call zonelines (zone5x , zone5y , ntop , nvs , NLINE , 80 ) 
Cc 
c 
c 

c c ********Zone 6 

c Write out Boundaries for Zone 6 
c 

write(80 , 96)6 ,nds+l , nvs , * ***Zone 6 ***’ 

c Zone 6 Line 1 

NLINE=1 

DATA (IC(6,I,1), 1 = 1 , 8 )/ 1,0, 5, 3, 4, 1,0,0/ 
write(80 , 97 ) (IC( 6 , I, NLINE) ,1=1,8) 
c 
C 

c Zone 6 Line 2 

NLINE=2 

DATA (IC(6,I,2),I=l,8)/2,6,0,0,6,2,0,0/ 
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c 



c 

c 



c 

c 



c 



c 



write(80 ,97 ) (IC(6, I,NLINE) , 1=1,8) 

write(80,*)zone5x(nbot-l ,1) ,zone5y(nbot-l , 1) 
call zonelines(zone6x,zone6y ,nds , nvs, NLINE, 80) 

Zone 6 Line 3 
NLINE=3 

DATA (IC(6,I,3),I=l,8)/3,2,0,0,6,3,0,0/ 
write(80 , 97 ) (IC(6 , I ,NLINE) ,1=1,8) 

call zonel ines ( zone6x , zone6y , nds , nvs , NLINE , 80) 

Zone 6 Line 4 

NLINE=4 

DATA (IC(6,I,4) ,1=1,8)74,1,0,0,6,4,0,0/ 
write(80 , 97) (IC(6 , I ,NLINE) , 1=1,8) 

write (80,*) zone5x(nbot-l ,nvs) ,zone5y (nbot-1 ,nvs) 
call zonel ines (zone6x , zone6y , nds , nvs , NLINE , 80) 

CL0SE(80) 



IN ORDER TO WRITE THE RESTART FILE, ZONES 2, 3, 5, 6 NEED TO 
BE APPENDED WITH THE "N-l" COLUMN FROM THE UPSTREAM ARRAY 
SO THAT OVERLAPPING CONTINUITY CAN BE MAINTAINED 
Append zone 2 

CALL MAKEDUM (ZONE IX , ZONE 1 Y , NUS , NVS , Z0NE2X , Z0NE2Y , NTOP , NVS , 
+DZ0NE2X , DZ0NE2Y , NTOPD , NVSD ) 

Append zone 3 

CALL MAKEDUM (Z0NE2X , Z0NE2Y , NTOP , NVS , Z0NE3X , Z0NE3Y , NDS , NVS , 
+DZ0NE3X , DZ0NE3Y ,NDSD ,NVSD) 

Append zone 5 

CALL MAKEDUM ( Z0NE4X , Z0NE4Y , NUS , NVS , Z0NE5X , Z0NE5Y , NBOT , NVS , 
+DZ0NE5X , DZ0NE5Y , NBOTD , NVSD) 

Append zone 6 

CALL MAKEDUM (ZONESX , Z0NE5Y , NBOT , NVS , Z0NE6X , Z0NE6Y , NDS , NVS , 
+DZ0NE6X , DZ0NE6Y , NDSD ,NVSD ) 



******************************************************************** 
WRITE INMESH RESTART FILE 
IF (NUMIT.GT. 1)THEN 

WRITE (* , ’ ( A) ’ ) 'WRITING FILE: ’ // RESTFILE 

OPEN (UNIT=78, FILE=RESTFILE,FORM=’ UNFORMATTED ’ , 

+ STATUS=’ UNKNOWN’) 

NZ0NE=6 

WRITE (78 )NZ0NE 
****Z0NE 1******** 

NZ0NE=1 

WRITE(78)NUS ,NVS 

WRITE(78) ((Z0NE1X(I , J) ,Z0NE1Y(I,J) ,1=1, NUS) , J=1,NVS) 

WRITE (78) ((IC(NZONE,M,N) ,M=1,8) ,N=1,4) 

****Z0NE 2******** 

NZ0NE=2 

WRITE (78) NTOPD, NVSD 

WRITE(78) ((DZ0NE2X(I, J) ,DZ0NE2Y(I , J) , 1=1, NTOPD), J=l, NVSD) 
WRITE(78) ((IC(NZ0NE,M,N),M=1,8),N=1,4) 

****Z0NE 3******** 

NZ0NE=3 

WRITE (78 )NDSD, NVSD 

WR ITE ( 78 ) ( (DZ0NE3X ( I , J ) , DZ0NE3Y ( I , J ) , 1= 1 , NDSD ) , J= 1 , NVSD ) 
WRITE(78) ( (IC(NZONE,M,N) ,M=1 ,8) ,N=1 ,4) 
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C ****ZQNE 4******** 

NZ0NE=4 

WRITE(78)NUS ,NVS 

WRITE ( 78 ) ( ( Z0NE4X (I , J) , Z0NE4Y ( I , J ) , 1= 1 , NUS ) , J=1 , NVS ) 
WRITEC78) ((IC(NZ0NE,M,N),M=1,8),N=1,4) 

C 

C ****ZONE 5******** 

NZ0NE=5 

WRITE(78)NB0TD ,NVSD 

WRITE(78) ( (DZ0NE5X(I , J) .DZ0NE5Y (I , J) ,1=1 ,NBOTD) , J=1,NVSD) 
WRITE (78) ((IC(NZONE,M,N) ,M=1,8) ,N=1,4) 

C 

C ****ZONE 6******** 

NZ0NE=6 

WRITE(78)NDSD,NVSD 

WRITE(78) ( (DZ0NE6X(I , J) ,DZ0NE6Y(I , J) , 1=1 ,NDSD) , J=1 ,NVSD) 
WRITE(78) ( (IC(NZONE,M ,N) ,M=1 , 8) ,N=1 ,4) 

C 

CL0SE(78) 

ENDIF 

ELSE 

write(* , *) * INMESH input files have not been written.* 

ENDIF 



C*** ******************************************************** ************ 
WRITE(* , *) *FIT2D COMPLETE* 

END 

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 

c 

C This subroutine takes rans wake data and makes it in the defining 
c line for the downstream grid zone boundary, 
c 

SUBROUTINE RANSWAKE(DSL , XTE , YTE, 

+ XYIN, XOUT , YOUT , NOUT ,NUMIN ) 

IMPLICIT NONE 

REAL DSL, XTE, YTE, PI , TEMPI ,TEMP2 

INTEGER NOUT,NUMIN,NI,I,K 

PARAMETER (NI=200 ,PI=3 . 141592653589793E00) 



C 

C ARRAYS 

C 

REAL XYIN(2,NI) , X0UT(NI ) , Y0UT(NI ) , YTEMP(NI) 

REAL XTEMP(NI) ,WAKCUB(4*NI) 
c 

C DEBUGGER WHAT IS COMING IN? 

C 

C WRITE (* , * ) (XYIN( 1 ,K) ,K=1 ,NUMIN) 

C WRITE (* , * ) ***** 

C WRITE(*,*)(XYIN(2,K),K=1,NUMIN) 

c CHECK THE DATA: 

C Does the stream line start aft of the trailing edge? 
c 

IF (XYIN(l , 1) .LT.XTE)THEN 

write (* , *) * ************************************ 1 
WRITE(* ,*)* ERROR: RANS wake starts forward of TE* 
write ( * , *) * ************************************ * 
NUMIN=-1 
ENDIF 
c 

c Does the stream data end before the domain runs out? 
TEMP1=XTE+DSL 

IF (XYIN(1,NUMIN).GT. TEMPI) THEN 
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wr it e(* , *) ’ ************************************ » 

WRITE (*,*) 'ERROR: RANS wake extends past flow domain 1 
write (* , *) ’ ************************* ************ 
NUMIN=-1 
ENDIF 
c 

C IF THE DATA IS OK THEN PROCEED THROUGH THIS: 

IF (NUMIN.GT.O)THEN 
C 

C FIRST THE LINE NEED TO BE FIXE UP A LITTLE. 

C THE BEGINNING NEEDS TO START AT THE TRAILING EDGE 

C OF THE FOIL AND THE END NEED TO BE AT THE DOWN 

C STREAM LIMIT OF THE FOIL 
C 

DO 100, 1=2 ,NUMIN+1 
XTEMP(I)=XYIN(l ,1-1) 

YTEMP ( I ) =X YIN (2,1-1) 

100 CONTINUE 
C 

C PUT BEGINNING AT THE TRAILING EDGE 
XTEMP( 1)=XTE 
YTEMP (1)=YTE 
C 

C PUT THE END AT THE DOWNSTREAM LIMIT 
XTEMP (NUMIN+2)=TEMP1 
YTEMP (NUMIN+2 ) = YTEMP (NUMIN+ 1 ) 

C 

C SPLINE THE LINE AND GET THE CUBIC COEFFICIENTS 
CALL UGLYDK (NUMIN+2 ,1,1, XTEMP , YTEMP ,0,0 , WAKCUB ) 

C 

C DEBUGGER 

C WRITE(* , *) ’ ONE OF THE CUBICS IS : ’ , WAKCUB (5) 

C WE WANT TO SEND NOUT X,Y POINTS BACK TO MAIN 
C DSL IS HOW FAR TO MARCH TOTAL 

TEMP 2=D SL/ ( FLOAT ( NOUT ) - 1 . 0 ) 

C DEBUGGER 

C WRITE (*,*) ’INTERVAL IS: , ,TEMP2 

X0UT(1)=XTE 

DO 200 1=2, NOUT 

XOUT ( I ) =XOUT (1-1) +TEMP 2 
C DEBUGGER 

C WRITE(* , *) ’THE XOUT VALUE IS:’,XOUT(I) 

200 CONTINUE 
C 

C NOW GET THE CORRESPONDING Y VALUES 

CALL EVALDK ( NUM IN+2 , NOUT , XTEMP , XOUT , YOUT , WAKCUB ) 

C 

C 

WRITE (* , *) ’ *** ’ 

WRITE(*,*)’***RANS data incorporated IN TO GRID***’ 
WRITE(*,*) ’***’ 
c WRITE ( * , * ) NUMIN , NOUT 

c WRITE(* , * ) ( (XOUT(k) , YOUT(k) ) ,K=1,N0UT) 

ENDIF 



RETURN 

END 

C 

C 

c 

C This subroutine uses a two-d vortex lattice lifting line to 
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c determine find the wake dividing line behind an arbitrary 

c foil section. This subroutine is a generalization of the 

c VLM2D code presented in the 13.04 course notes. 

C 

SUBROUTINE ADAPT (PRESSX , PRESSY , SUCTX , SUCTY , WAKEX , WAKEY , 

+ NPRES , NSUCT , DSL , NWAKE , TUNPARAM ) 

IMPLICIT NONE 

INTEGER NPRES , NSUCT , NWAKE ,I,J,K,N,M,NI, NPTS , NP AN , IERR , NTQT 
INTEGER NIMAGEPR 
REAL PI 

PARAMETER (NI=200 , PI=3 . 14 1592653589793E00) 

REAL DSL , DELS , LEN , DX , DY , TOP , s urn , ENDX , U , V , TEMP , STEP 
REAL TUNPARAM, RTEMP , DXR , DYR, YIMAGE, signimage 
C 

C ARRAYS 
C 

REAL PRESSX (NI ) , PRESSY (NI ) , SUCTX (NI ) , SUCTY (NI ) 

REAL WAKEX(NI) , WAKEY (NI ), cambx(ni) ,camby(ni) 

REAL T0PCUBX(4*NI) ,B0TCUBX(4*NI) , CAMBCUBX(4*NI) 

REAL TOPCUBY (4*NI) ,B0TCUBY(4*NI) , CAMBCUBY(4*NI ) 

REAL WAKCUB (4*NI) 

REAL ARCTOP (NI ) , ARCBOT (NI ) , STEMP (NI ) , ARCAMB (NI ) 

REAL TEMPTX(NI) ,TEMPTY(NI) 

REAL TEMPBX (NI ) , TEMPBY (NI ) 

REAL SV(NI) ,SC(NI) ,DS(NI) 

REAL XV (NI ) , YV (NI) ,XC(NI) , YC(NI) , B(NI) ,TX(NI) ,TY(NI) 

REAL XN(NI) , YN(NI) ,A(NI ,NI) ,GAMMA(NI) ,WKAREA(NI) ,dydx(ni) 

INTEGER IPIVOT(NI) 
c 
c 

c What is the Y distance to the first image set? ) 

c (the "Image Plane" is at half this distance 

YIMAGE=1/ (TUNPARAM) 
c 

c write(* ,*) ’the wall is at : 9 , ywall 

C Spline Foil Offsets and determine cubic coefficients, 

c 

C Array ARC is returned with arclength params non-dimed from 0..1 

C 

CALL PUGLYDK (NSUCT , SUCTX , SUCTY , ARCTOP , TOPCUBX , TOPCUBY ) 

CALL PUGLYDK (NPRES , PRESSX , PRESSY , ARCBOT , BOTCUBX , BOTCUBY ) 
c 
c 

c Extract Points along the top and bottom surface to find the 

c mean camber line, 

c 

C HOW MANY POINTS TO EXTRACT? 

NPTS=40 

C 

c 

C AT WHAT ARC COORDINATE LOCATIONS TO EXTRACT? 

C 

c write(*,*) arctop(nsuct ) 

STEMP (1)=0.0 
DO 100 1=1, NPTS 

STEMP ( I ) = (FLOAT ( I- 1 ) /FLOAT (NPTS-1 ) )* ARCTOP (NSUCT ) 
c write(* ,*)stemp(i) 

100 CONTINUE 
C 

C EXTRACT THE VALUES ON THE SUCTION SIDE: 

C 

CALL PEVALDK (NSUCT , NPTS , STEMP , ARCTOP , TEMPTX , TEMPTY , TOPCUBX 
+, TOPCUBY) 
c 
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c 

C AT WHAT ARC COORDINATE LOCATIONS? 

STEMP( 1)=0 . 0 

DO 200 1=1 ,NPTS 

STEMP ( I ) = ( FLO AT ( I - 1 ) / FLO AT ( NPTS - 1 ) ) * ARCBOT ( NPRES ) 

200 CONTINUE 
C 

C EXTRACT THE VALUES ON THE PRESSURE SIDE: 

C 

CALL PEVALDK ( NPRES , NPTS , STEMP , ARCBOT , TEMPBX , TEMPBY , BOTCUBX 
+,BOTCUBY) 
c 

C FIND THE MEAN CAMBER LINE(SIMPLE MEAN OF TWO COORDINATES): 
c 

DO 300 1=1, NPTS 

C AMBX ( I ) = (TEMPTX ( I )+TEMPBX (NPTS-I+1 ) ) *0 . 50 
C AMBY ( I ) = (TEMPTY ( I ) +TEMPBY (NPTS-I+1 ) ) *0 . 50 
C 

C debugger 

C WRITE (88 , *) CAMBX (I ) , CAMBY(I ) 

300 CONTINUE 
c debugger 

c 

c Spline the mean camber line to obtain the Cubic Coefficients 

c Once again using the arc parameter scaling. 

CALL PUGLYDK (NPTS , CAMBX , CAMBY , ARCAMB , CAMBCUBX , CAMBCUBY) 

C 

C DETERMINE THE LOCATIONS OF THE VORTICES AND CONTROL POINTS 
C 

C How many vortices and control points should there be? 
c A convergence study was performed by varying the number 

c of panels. The difference in lift coefficient is less 

c than 0.1*/, when stepping from 20 to 40 panels. So. . . 

c 20 panels is sufficient. 

NPAN=40 

c 

c How many pairs of image foils should there be? 
c Obviously, the images should only be put in in pairs, 

c This value will be adjusted until there is a converged 

c lift coefficient, 

c 

c A convergence study was conducted. The relative error is 

c approx 0.04*/, with 16 image pairs. This is sufficient 

c for this application. 

NIMAGEPR=16 

c 

c Establish COSINE spacing on ARC Length Parameter 

c this is straight out of the 13.04 notes, 

c DELS is the COSINE spacing interval 

DELS=PI/NPAN 

DO 400 1=1 ,NPAN 

S V ( I ) =0 . 50 * ( 1 . 0-COS ((1-0.50) *DELS ) ) 

SC (I )=0 . 50* ( 1 . 0-C0S(I*DELS) ) 

DS( I )=PI*SQRT(SV(I )* ( 1 .O-SV(I) ) )/NPAN 
400 CONTINUE 
C 

C EXTRACT THE X AND Y VALUES OF THE VORTICES AND CONTROL POINTS 
C 

CALL PEVALDK (NPTS , NP AN , S V , ARCAMB , XV , YV , CAMBCUBX 
+, CAMBCUBY) 

CALL PEVALDK (NPTS , NP AN , SC , ARC AMB , XC , YC , CAMBCUBX 
+, CAMBCUBY) 
c 
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To Panel the walls insert that stuff here: 



c 
c 

c How far away are the walls: TUNPARAM dictates 

c How far do you want to panel up an down stream of 

c the foil? DSL is already passed in. So a reasonalbe 

c limit for paneling is LE - DSL to TE + DSL + 1 

c 

c What is a reasonable spacing: UNIFORM 

c How many panels per WALL: 2*NPAN seems OK 

c 

c write locations out to datafile for inspection to fort. 88 
c debugger 

WRITE ( 88 , * ) ' VARIABLES=XV , YV , XC , YC ' 

WRITE(88 , *) 'ZONE T=V0RT' 

DO 500 1=1 ,NPAN 

WRITE(88, * )XV(I) ,YV(I),XC(I),YC(I) 

500 CONTINUE 

C 

c Calculate the UNIT NORMAL XN,YN at each control point 
c unit normal not used for anything yet, but it may be 
c useful in the future, 

c 

c Also, slope is required on the foil for V*n equation later 
c 

c write(88 , *) 'zone ' 

DO 600 1=1, NP AN-1 
DX=XV ( 1+1 ) -XV ( I ) 

DY=YV(I+1)-YV(I) 

LEN=SQRT(DX*DX+DY*DY) 

XN ( I ) =-DY/LEN 

YN(I)=DX/LEN 

dydx(i)=dy/dx 

c write(88,*) xn(i),yn(i) 

600 CONTINUE 

FIX UP THE VALUE AT THE TRAILING EDGE BECAUSE THERE IS NO 
VORTEX AFTER THE LAST CONTROL POINT. 

DX=XC(NPAN)-XV (NPAN) 

DY=YC (NPAN ) - YV (NPAN ) 

XN(NPAN)=-DY/SQRT (DX*DX+DY*DY ) 
YN(NPAN)=DX/SQRT(DX*DX+DY*DY) 
dydx ( npan ) =dy / dx 
c write(88, * )xn(npan) ,yn(npan) 

C 

c ALL THE UNIT NORMALS & Slopes ARE NOW ESTABLISHED. 

C 

C COMPUTE THE INFLUENCE COEFFICIENTS A(N,M) AND THEN 
C INVERT THE MATRIX . (linearized) 

C The A matrix are the influence coefficients. Based on 
c poor agreement with RANS, it was decided to go to the 

c true locations of the vortices on the mean camber line 

c vice a linearized surface. 



T0P=1 . 0/(2 . 0*PI) 

DO 700 N=1 ,NPAN 
DO 700 M=1 ,NPAN 

rtemp=sqrt ( (xv(ra)-xc(n) ) **2+(yv(m)-yc(n) )**2) 

dyr=-(xc(n)-xv(m))/rtemp 

dxr=-(yc(n)-yv(m) )/rtemp 
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A(N,M)=(TOP/(rtemp) ) * (dxr*xn(n)+dyr*yn(n) ) 
c write (89 , *)a(n ,m) 

c The above include the influence of the vortex influence on 

c the foil itself. Now, add in the images above and below 

c the foils due to wall effects. Need Tunparam 
c 

DO 710 K=1 ,NIMAGEPR 
c 

c Test K to see if it is a positive or negative image line 
c if true then it is a positive image line 

c if not true then it is a negative image 

c 

if ( (2*int (K/2) ) .eq.k) then 
signimage=l .0 
else 

signimage=-l . 0 
endif 

c DEBUGGER 

c if ( (n . eq . 1 ) . and . (m . eq . 1 ) ) then 

c write(*,*) signimage 

c endif 

c First Add in the Effect due to the image ABOVE (in pos y-dir) 
c 

rtemp=sqrt( (xv(m)-xc(n) )**2+ 

'/. ( (k*yimage+ signimage*yv(m) )-yc(n) )**2) 

dyr=-(xc(n)-xv(m) )/rtemp 

dxr=-(yc (n)-(k*yimage+signimage*yv(m) ) )/rtemp 
A(N,M)= A(N,M)+signimage*(TOP/(rtemp) ) * (dxr*xn(n)+dyr*yn(n) ) 
c 

c Next, Add in the Effect due to the image BELOW (in neg y-dir) 
c 

rtemp=sqrt( (xv(m)-xc(n) )**2+ 

'/. ( (-k*yimage+ signimage*yv(m) )-yc(n) )**2) 

dyr=-(xc(n)-xv(m) )/rtemp 

dxr=-(yc(n)-(-k*yimage+signimage*yv(m) ) )/rtemp 
A(N,M)= A(N,M)+signimage*(TOP/(rtemp) )*(dxr*xn(n)+dyr*yn(n) ) 

710 continue 

c 

c 

c 

700 CONTINUE 

FACTOR THE A MATRIX 

CALL FACT0R( A , IPI V0T , WKAREA , NPAN , NI , I ERR) 

COMPUTE THE w VELOCITY AT THE NTH CONTROL POINT 
THE B MATRIX OF A*GAMMA=B 

DO 800 N=1 ,NPAN 
B(N)=DYDX(N) 
write (* , *)dydx(n) 

CONTINUE 

NOW BACK SUBTITUTE TO EXTRACT THE VORTEX STRENGTHS 

CALL SUBST( A , B , GAMMA , IPIV0T , NPAN , NI ) 

check the circulation distribution on the foil 
c 

c write out to fort. 89 

c 
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sum=0 . 0 

do 900 i=l,npan 
s um= s um+g amma ( i ) 
write(89 , *)xv(i) ,gamma(i)/ds (i) 

900 continue 

write(*,*) ’Estimated CL=’,sum*2.0 
write(*,*) ’Using Vortex Lattice’ 
c 

c extract the wake 

c 

c The circulation is now known for each vortex point, 
c To find the wake, "Step 11 off the trailing edge and 

c evaluate the field point velocity. March off a distance 

c in that direction. Evaluate the field point velocity again, 

c adjust course. Keep a record of the path. That will be 

c the line for the wake centerline, 

c 

c The following quantities are needed: 

c NWAKE, DSL 

c Want to calculate: WAKEX, WAKEY <Q NWAKE points between TE and DSL 

c 

c March down the wake first, spline the values and then extract 
c the NWAKE points, 

c 

c Where is the Trailing edge? <0 pressx(l) ,pressy (y) 

c 

c Where do you want to stop? 0 pressx(l)+DSL & corresponding y 
c 

endx=pressx(l)+dsl+l .0 

WHAT STEP DO YOU WANT TO MARCH IN? (think of as scaled time) 
STEP=DSL/25 . 0 

START MARCHING 

1=1 

TX ( 1 ) =PRESSX ( 1 ) 

TY ( 1 ) =PRESSY ( 1 ) 
write(88,*) ’zone’ 

DO 910 WHILE (TX(I) .LT.ENDX) 

FIND FIELD POINT VELOCITY DISTURBANCE DUE TO VORTICES 
REFERENCE NEWMAN PAGE 190 
c U=1 to add in the freestream effect 
U=1.0 
V=0.0 

DO 920 N=1 ,NPAN 

TEMP=2*PI* ( (XV(N)-TX(I ) )*(XV(N)-TX(I) )+(YV(N)-TY(I ) )* 

V. (YV (N)-TY(I) ) ) 

U=U+GAMMA (N) *(TY(I)-YV(N) ) /TEMP 
V=V-G AMMA (N) * (TX ( I ) -XV (N ) ) /TEMP 
c 

c Add in the influence of the images 

c 

c 

DO 1950 K=1 ,NIMAGEPR 
c 

c Test K to see if it is a positive or negative image line 
c if true then it is a positive image line 

c if not true then it is a negative image 

c 

if ( (2*int (K/2) ) .eq.k) then 
signimage=l . 0 
else 

signimage=-l .0 
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endif 
c DEBUGGER 

c First Add in the Effect due to the image ABOVE (in pos y-dir) 
c 



temp=2*PI*( ( (XV(N)-TX(I) )**2+ 

’/• ( (k*yimage+ signimage*yv(n) )-ty(i) )**2) ) 

U=U+s ignimage*GAMMA (N) * (TY ( I ) - (k*yimage+ 

*/. signimage*yv(n) ) )/TEMP 

V=V-s ignimage*GAMMA (N) * (TX ( I ) -XV (N) ) /TEMP 
c 

c Next, Add in the Effect due to the image BELOW (in neg y-dir) 



temp=2*PI*( ( (XV(N)-TX(I) )**2+ 

*/• ( (-k*yimage+ s ignimage*yv(n) )-ty (i) ) **2) ) 

U=U+signimage*GAMMA(N)* (TY(I)-(-k*yimage+ 

% signimage*yv(n) ) )/TEMP 

V=V-s ignimage+GAMMA (N) * (TX ( I ) -XV (N) ) /TEMP 



1950 continue 



c 

c 

c 

c 

c 

920 continue 

CALCULATE THE NEXT LOCATION TO LOOK 

Make the first step a baby step 

IF(I -LE. 3) THEN 
STEP=STEP/8 . 0 
ELSE 

STEP=DSL/25 . 0 
ENDIF 



Where is the next uake location? 

TX ( 1+ 1 ) =TX ( I ) +U*STEP 
TY ( 1+1 ) =TY ( I ) +V+STEP 
write (88 , *)tx(i) ,ty(i) ,0,0 
1 = 1+1 

910 END DO 
NT0T=I 

SPLINE THE POINTS TO OBTAIN CUBIC COEFFICIENTS 

...SPLINE SPACING WITH FIXED SLOPE AT THE ENDS 
CALL UGLYDK ( NTOT , 1 , 1 , TX , TY , 0 , 0 , WAKCUB ) 

...EVALUATE SPLINE TO FIND STEP SIZE AT INTERMEDIATE POINTS 

WHERE ARE THE OUTPUT POINTS? 

WAKEX ( 1 ) =PRESSX ( 1 ) 

DO 950 I=1,NWAKE-1 

WAKEX ( 1+ 1 ) =WAKEX ( I ) +D SL/FLOAT ( ( NWAKE- 1 ) ) 

950 CONTINUE 

EXTRACT THE POINTS FOR THE WAKEX AND WAKEY TO RETURN TO MAIN 
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CALL EV ALDK (NTOT , NWAKE ,TX , WAKEX , WAKEY , WAKCUB ) 
write (88,*) 'ZONE T=WAKEADS ' 

DO 960 1=1,10 

WRITE (88 , *) WAKEX(I) ,WAKEY(I) ,0,0 
960 CONTINUE 

write (* ,*) ' ************************************************** 

write(* , *) 'THE GRID BOUNDARIES HAVE BEEN ADAPTED TO THE WAKE' 

write(*,*) » USING VORTEX LATTICE METHOD' 

write (* ,*) >************************************* ************* 

RETURN 

END 



************************************************************************ 
* SINGLE PRECISION VERSION OF DAVE GREELEY'S DIRECT SOLVER * 

************************************************************************ 
C234567890 1234567890 1234567890 1234567890 12345 67 890 1234567890 1234567890 12 

C Subroutine FACTOR 

SUBROUTINE FACTOR(W, IP IVOT , D , N , NDIM , IERR) 

IMPLICIT REAL(A-H,0-Z) 

DIMENSION W(NDIM,NDIM) ,IPIVOT(*) ,D(*) 

IERR=0 
DO 10 1=1, N 
IPIVOT(I)=I 
R0WMAX=0 . 

DO 9 J=1,N 

ROWMAX=MAX (ROWMAX , ABS ( W ( I , J ) ) ) 

9 CONTINUE 

IF (ROWMAX. EQ. 0. ) GO TO 999 
D ( I ) =ROWMAX 

10 CONTINUE 
NM1=N-1 

IF(NMl.EQ.O) RETURN 
DO 20 K= 1 , NM 1 
J=K 

KP1=K+1 
IP=IPIVOT (K) 

COLMAX=ABS ( W ( IP , K ) ) /D ( IP ) 

DO 11 I=KP1,N 
IP=IPIVOT (I) 

AWIK0V=ABS(W(IP,K))/D(IP) 

IF(AWIKOV.LE.COLMAX) GO TO 11 

COLMAX=AWIKOV 

J=I 

11 CONTINUE 

IF (COLMAX . EQ . 0 . ) GO TO 999 
IPK=IPIVOT(J) 

IPIVOT(J)=IPIVOT(K) 

IPIVOT (K) =IPK 
DO 20 I=KP1,N 
IP=IPIVOT(I) 

W(IP ,K)=W(IP ,K)/W(IPK,K) 

RATIO=-W(IP,K) 

DO 20 J=KP1,N 

W(IP, J)=RATIO*W(IPK, J)+W(IP , J) 

20 CONTINUE 

IF(W(IP,N) .EQ.O.) GO TO 999 
RETURN 
999 IERR=2 
RETURN 
END 

C End of FACTOR 
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C Subroutine SUBST 

c* ******************************************************* *************** 



SUBROUTINE SUBST(V,B , X , IPIVOT , N ,NDIM) 

IMPLICIT REAL(A-H,0-Z) 

DIMENSION W(NDIM,NDIM) ,B(*) ,X(*) , IPIV0T(*) 

IF(N.GT.l) GO TO 10 
X(l)=B(l)/W(l,l) 

RETURN 

10 IP=IPIV0T(1) 

X(1)=B (IP) 

DO 15 K=2,N 
IP=IPIV0T(K) 

KM1=K-1 
SUM=0 . 

DO 14 J=1 ,KM1 

SUM=W(IP , J) *X(J)+SUM 

14 CONTINUE 
X(K)=B(IP)-SUM 

15 CONTINUE 
X(N)=X(N)/W(IP,N) 

K=N 

DO 20 NP1MK=2,N 
KP1=K 
K=K-1 

IP=IPIVOT(K) 

SUM=0 . 0 
DO 19 J=KP1,N 

SUM=W(IP , J)*X(J)+SUM 

19 CONTINUE 
X(K)=(X(K)-SUM)/W(IP,K) 

20 CONTINUE 
RETURN 
END 

C End of SUBST 

************************************************************************ 
♦END OF SINGLE PRECISION VERSION OF DAVE GREELEY’S DIRECT SOLVER * 

c 

C+++++++++++++++^+++-^+++-H-++++++++-H-+++-H-++-H-+++-H-++++-f+++-H-+++-H-++++4*++ 

C 

SUBROUTINE MAKEDUM (AX , AY , IMAXA , JMAXA , BX , B Y , IMAXB , JMAXB , 

+X0UT , YOUT , IMAXOUT , JMAXOUT ) 



This subroutine takes in two axially adjacent arrays and 
appends the imaxa-1 vertical line of points to the left 
hand side of the downstream array. 

Last Updated: 11 November 

IMPLICIT NONE 

INTEGER IMAXA , JMAXA , IMAXB , JMAXB , IMAXOUT , JMAXOUT 
INTEGER I ,NI , J 
PARAMETER (NI=200) 

ARRAYS 

REAL AX(NI.NI) ,AY(NI,NI) ,BX(NI ,NI) ,BY(NI,NI) 

REAL X0UT(NI,NI) ,Y0UT(NI,NI) 
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DO 100 J=1 , JMAXA 

XOUT ( 1 , J)=AX(IMAXA-1 , J ) 

YOUT ( 1 , J ) =A Y ( IMAXA-1 , J ) 

100 CONTINUE 

DO 200 1=1 , IMAXB 
DO 300 J=1 , JMAXB 
X0UT(I+1 , J)=BX(I , J) 

Y0UT(I+1,J)=BY(I,J) 

300 CONTINUE 
200 CONTINUE 

IMAX0UT=IMAXB+1 
JMAXOUT= JMAXB 
RETURN 
END 

C+++++++++++++++4^+++^+++^++++++++-M'+++-M'++-Hf+++-M-+++-J-++++-M-+++-M'+++++++ 

SUBROUTINE ZONELINES ( X , Y , IMAX , JMAX , NLINE , NFILE ) 

C 

C This subroutine assists in writing output file for FIT2D 
c to INMESH output file, 

c 

c X,Y are arrays containing points 

c IMAX, JMAX are # of array elements 
c NLINE is the line number in the zone 

c NFIEL is the file output number 

IMPLICIT NONE 

INTEGER NI ,NJ , I , IMAX, JMAX , J , NLINE, NFILE 
PARAMETER (NI=200 ,NJ=200) 

REAL X(NI,NJ),Y(NI,NJ) 

REAL X 1 (NI , N J ) , X2 (NI , N J ) , Y 1 ( NI , N J ) , Y2 (NI , N J ) 

REAL RI 1 ,RI2 ,RJ 1 ,RJ2 ,BETA 
C 
C 
C 

IF (NLINE. EQ.l) THEN 
DO 100 J=1 , JMAX 

WRITE (NFILE, *)X( 1,J) ,Y(1,J) 

100 CONTINUE 
ENDIF 

IF (NLINE. EQ. 3) THEN 
DO 300 J=1 , JMAX 

WRITE (NFILE, *)X( IMAX , J) ,Y(IMAX,J) 

300 CONTINUE 
ENDIF 

IF (NLINE. EQ. 2) THEN 
DO 200 1=1, IMAX 

WRITE(NFILE, *)X(I , l) ,Y(I,l) 

200 CONTINUE 
ENDIF 

IF (NLINE. EQ .4) THEN 
DO 400 1=1, IMAX 

WRITE(NFILE,*)X(I, JMAX), Y(I, JMAX) 

400 CONTINUE 
ENDIF 



RETURN 

END 



C 

SUBROUTINE B0RDERS(X,Y,X1,Y1,X2,Y2,X3,Y3,X4,Y4,IX,JY,NI) 

C 

C This subroutine takes in the x,y points which represent the four 
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c 

c 
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c 

c 

c 

c 



lines surrounding a 2d box. It then writes them to a 2d array 
which represents the box boundaries. The "interior" of the 2d 
array will still be zeros. The interior is filled in another 
subroutine which interpolates between all the points. 



Side 1 



Side 4 
ZONE "n" 

Side 2 



Side 3 



/\ 

I Array reference to enter finite 

I interpolation scheme. 

JY I 

I > 

IX 



IMPLICIT NONE 

INTEGER I, IX, JY,NTEMP,NI 

ARRAYS 

REAL X(NI,NI) , Y(NI ,NI) ,X1(NI) , Y1 (NI) ,X2(NI) , Y2(NI) 

REAL X3(NI) ,Y3(NI) ,X4(NI) ,Y4(NI) 

REAL X5(NI) ,Y5(NI) ,X6(NI) ,Y6(NI) 

Check line One: Determine if the order forward up or forward down. 

ntemp = 0 then count forwaxd 
ntemp = 1 count backwards 

debugger 

write(* , *) ’Debugging in borders y(jmax), y(l)’ 
write(* , *) yl(jy),yl(l) 

NTEMP=0 

IF (Yl(JY) .LT. Yl(l) )THEN 
NTEMP=1 
ENDIF 

Now put line one in to x and y arrays 



DO 100 1=1, JY 

IF (NTEMP. EQ.O) THEN 
X(1,I)=X1(I) 

Y(1 , I)=Y1 (I) 

ELSE 

X( 1 , I )=X1 ( JY-I+1 ) 
Y(1 , 1)=Y 1 (JY-I+1 ) 
ENDIF 

debugger 



if (i.eq.l)then 

write(80, *) "ZONE" 
endif 

write(80,*) x(l,I) ,y(l,I) 
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c end debugger 
c 

100 CONTINUE 
c 

Check line Two: Determine if the order is forward or backwards 

ntemp = 0 then count forward 
ntemp = 1 count backwards 
NTEMP=0 

IF (X2(IX) .LT.X2(1))THEN 
NTEMP=1 
ENDIF 

Now put line two in to x and y arrays 

DO 200 1=1, IX 

IF (NTEMP. EQ.O) THEN 
X(I,1)=X2(I) 

Y(I , l)=Y2(I) 

ELSE 

X(I,1)=X2(IX-I+1) 

Y(I,1)=Y2(IX-I+1) 

ENDIF 

debugger 

if (i.eq.l)then 
write(80, *) “ZONE" 
endif 

write(80, *) x(I , 1) ,y(I , 1 ) 
end debugger 
CONTINUE 

Check line Three : Determine if the order forward up or forward down. 

ntemp = 0 then count forward 
ntemp = 1 count backwards 

NTEMP=0 

IF (Y3(JY) .LT.Y3(1) )THEN 
NTEMP=1 
ENDIF 

Now put line three in to x and y arrays 

DO 300 1=1, JY 

IF (NTEMP. EQ.O) THEN 
X(ix,I)=X3(l) 

Y(ix,I)=Y3(I) 

ELSE 

X(ix, I)=X3( JY-I+1) 

Y(ix,I)=Y3(JY-I+l) 

ENDIF 

debugger 
c 

c if (i.eq.l)then 

c write(80,*) “ZONE" 

c endif 

c write(80,*) x(ix, I) ,y(ix , I) 

c 

c end debugger 

c 
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300 CONTINUE 
c 

C* * * * 



Check line Four: Determine if the order is forward or backwards 

ntemp = 0 then count forward 
ntemp = 1 count backwards 
NTEMP=0 

IF (X4(IX).LT.X4(1))THEN 
NTEMP=1 
END IF 

Now put line four in to x and y arrays 

DO 400 1=1, IX 

IF (NTEMP. EQ.O) THEN 
X(I , jy)=X4(I) 

Y(I,jy)=Y4(I) 

ELSE 

X(I , jy)=X4(IX-I+l) 

Y(I, jy)=Y4(IX-I+l) 

ENDIF 



debugger 

if (i.eq.l)then 

write (80,*) "ZONE" 
endif 

write(80 , *) x(I , jy) ,y(I , jy) 
end debugger 
CONTINUE 



CALL interp(X,Y, IX, JY) 

RETURN 

END 



SUBROUTINE interp(X , Y , IMAX , JMAX) 



Below is a subroutine for doing ISOPARAMETRIC interpolation 
for a single zone. It assumes that all four edges axe 
already split up to their desired spacings. 

This method is based on the simple isoparemetric method presented 
c in "FINITE ELEMENT PROCEDURES" by Bathe, 
c 
c 

c Written by John Dannecker for use in FIT2D.F 

c Last Modified: November 7, 1996 

c 

IMPLICIT NONE 

INTEGER NI ,NJ , I , IMAX, JMAX, J 
PARAMETER (NI=200 ,NJ=200) 

REAL X(NI,NJ),Y(NI,NJ) 

REAL XI (NI ,NJ) ,X2(NI ,NJ) ,Y1(NI ,NJ) , Y2(NI ,NJ) 

REAL RI1 ,RI2 ,RJ 1 ,RJ2,beta 
c 

c i=15 

c j=25 

c DO ALL THE X(I,J) 
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c 



DO 100 1=2 , IMAX-1 

RI1=(X(I , 1)-X(1 , 1) )/(X(IMAX, 1)-X(1,1)) 
ri2=(x(i , jmax)-x(l , jmax) )/ (x(imax, jmax )-x(l , jmax) ) 
c RI1=RI1+(X(I , JMAX) -X(l , jmax) )/ (X(IMAX , JMAX)-X( 1 , JMAX) ) 

C 

DO 200 J=2 , JMAX-1 

beta=(Y(l ,J)-Y(1,1))/ (Y( 1 , JMAX)-Y(1 , 1) ) 
c RJ1=RJ 1+(Y(IMAX, J )-Y(IMAX, 1) )/ (Y(IMAX, Jmax)-Y(IMAX, 1)) 

X(I, J)=x(l, j)+(x( imax, j)-x(l, j) )*((ril*( 1.0-beta) )+(ri2* (beta) ) ) 

200 CONTINUE 

100 CONTINUE 
C 

C DO ALL THE Y(I,J) 

C 

DO 300 j=2 , jMAX-1 

Rjl=(y(l > j)-y(l, l))/(y(l, jmax)-y(l,l)) 
Rj2=(y(imax,j)-y(iraax,l))/(y(IMAX,JMAX)-y(imax,l)) 

C 

DO 400 i=2 , iMAX-1 

beta=(x(i, l)-x(l, 1) )/(x(imax, l)-x(l , 1) ) 
c Ri2=Ri2+(x(i, jmax )-x(l, jmax))/ (x(IMAX, jmax )-x(l, jmax)) 

Y(I,J)=y(i, l)+(y(i, jmax)-y(i, 1) )*(rj 1*(1 . 0-beta)+rj2*beta) 

400 CONTINUE 

300 CONTINUE 
RETURN 
END 
c 
c 

C+++++++++++++++-f^+++4H-+++4H-+++-H-+++-H-+++-f++++++++-H-+++-H-+++-H-+++-H-+++-H-++ 

C 

SUBROUTINE FINDYPLUS (RESBL , RESWALL , REYNOLD , TUNPARAM , USL , DSL) 
c 

c This subroutine calculates the parameter y+ 
c using the method in Anderson ,Tannehill ,Pletcher 
c Computaional Fluid Mechanics and Heat Transfer, 1984. 
c 

C Y+ is checked against the user input values for cell height 

c on the foil and wall. If the user cell spacing exceeds 

c y+ criteria, then cell spacing can be changed. 

IMPLICIT NONE 

REAL RESBL , RESWALL , REYNOLD , TUNPARAM , USL , DSL , HEIGHT , NU 
REAL LF , LW , VEL, YPLUSF , REYWALL , YPLUSW , TEMP 
CHARACTER ANSWER* 3 
C 

C MIT WATERTUNNEL CROSS SECTION HEIGHT (FEET) 

C 

HEIGHT=20. 0/12.0 
C 

C KINEMATIC VISCOSITY (NU FT~2/SEC) 

C FRESH WATER 0 70 DEGREES F. (SOURCE PNA 1967) 

NU=1 .0519E-5 
C 
C 

C CALCULATE CHORD LENGTH 

c 

LF=TUNPARAM*HEIGHT 

c 

c Calculate wall length 

c 

LW=TUNP ARAM* HEIGHT* ( 1 . 0+USL+DSL) 



c 
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c CALCULATE FLOW SPEED IN TUNNER (FT/SEC) 

C 

V EL=NU * RE YN 0 LD / LF 
C 

C CALCULATE YPLUS FOR FOIL 

C 

YPLUSF=(3 . 0/ (REYNOLD** (-0 . 10) *SQRT(0 . 0227) *VEL/NU) )/LF 
C 

C CALCULATE WALL BASED REYNOLD NUMBER 

C 

REYWALL=VEL*LW/NU 

C 

C CALCULATE YPLUS AT WALL 

C 

YPLUSW=(3 . 0/(REYWALL** (-0 . 10) *SQRT(0. 0227) *VEL/NU) )/LF 
C 

C COMPARE USER INPUT CELL HEIGHT TO Y+: 

C 

c 

c debugger 

c 

99 FORMAT (A 12 ,F10.6,A6,F10.6) 

write(*, *) ' Cell Dim User Input Val Yplus(3)' 
write(*,99) 'On Foil: ',resbl,' ',yplusf 
write(*,99) 'On Wall: ’,reswall,' ',yplusw 
c 

C ON THE FOIL: 

TEMP=YPLUSF/RESBL 
IF (RESBL . GE . YPLUSF ) THEN 

WRITE(* , *) 'USER INPUT CELL HEIGHT ON FOIL IS GREATER THAN' 
WRITE(* , *) 'CALCULATED Y+(3) VALUE.' 

WRITE(* , ' (A , $) ' ) 'DO YOU WANT TO CHANGE TO CALC Y+(3)?<n>:' 

READ(* , ' (A) ' ) ANSWER 

IF ((ANSWER. EQ. 'Y') .OR. (ANSWER. EQ. 'y')) THEN 
RESBL=YPLUSF 

WRITE(*,*) 'CELL HEIGHT ON FOIL SET TO Y+(3) . ' 

ELSE 

WRITE(*,*) 'CELL HEIGHT ON FOIL IS UNCHANGED.' 

END IF 
ELSE 

WRITE(*,*) 'USER SPECIFIED CELL HEIGHT ON FOIL IS ADEQUATE' 

ENDIF 

C 

C AT THE WALL 

TEMP=YPLUSW/RESWALL 
IF (RESWALL . GE . YPLUSW) THEN 

WRITE(* , *) 'USER INPUT CELL HEIGHT AT WALL IS GREATER THAN' 

WRITE(* , *) 'CALCULATED Y+(3) VALUE.' 

WRITE(* , ' (A,$) ' ) 'DO YOU WANT TO CHANGE TO CALC Y+(3)?<n>: ' 

READ(* , ' (A) ' ) ANSWER 

IF ((ANSWER. EQ. 'Y') .OR. (ANSWER. EQ. 'y')) THEN 
RESWALL=YPLUSW 

WRITE(*,*) 'CELL HEIGHT AT WALL SET TO Y+(3) . ' 

ELSE 

WRITE(*,*) 'CELL HEIGHT AT WALL IS UNCHANGED.' 

ENDIF 

ELSE 

WRITE(*,*) 'USER SPECIFIED CELL HEIGHT AT WALL IS ADEQUATE' 

ENDIF 

C 

C 

RETURN 

END 

C 

C++^+++++++++++-h-+++-w-+++^+++++++++++++-m-++-h-+++-h-++++++++-h-++++++++++++ 
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100 



c 



c 

200 



SUBROUTINE TUND (XBEG , YBEG , XEND , YEND ,RESBL , RESWALL , XOUT , YOUT , 
+NP0INTS, PACKING) 

LAST MODIFIED: 29 OCTOBER 1996 

This subroutine generates the vertical spacing of points on 
lines above and below the foil. 

IMPLICIT NONE 

REAL PI ,TW0PI , ZERO , ONE , HALF 

INTEGER NI , I , NNEARBL , NNEARWALL , NREMAIN 

PARAMETER (PI=3 . 14159 ,TV0PI=6 . 28318 , NI=200, ZER0=0.0, 0NE=1.0) 
INTEGER NPOINTS 

REAL XBEG , XEND , YBEG , YEND , RESBL , RESWALL , PACKING 
REAL DELX , DELY , TEMPX , TEMP Y , DLEFT , DRIGHT 
real xvect ,yvect ,hypo 

ARRAYS 

REAL XOUT (NI ), YOUT (NI) ,XLIN(NI) , YLIN(NI) , XREM(NI) ,YREM(NI) 



Determine spacing of points for boundary layer along foil 

How many points close packed near foil and at the wall? 

NNEARBL=INT(PACKING*FLOAT(NPOINTS) ) 

NNEARWALL=NNEARBL 

NREMAIN=NP0INTS-NNEARWALL-NNEARBL+2 



Size the cells on the foil: 

X0UT( 1)=XBEG 
Y0UT(1)=YBEG 
DELX=XEND-XBEG 
DELY=YEND-YBEG 

hypo=sqrt (delx*delx+dely*dely ) 

xvect=delx/hypo 

yvect=dely/hypo 

debugger 

write(* ,*)xvect ,yvect ,delx ,xend,xbeg 
Each cell is 25’/, larger than the previous cell. 
DO 100 1=1 , NNEARBL- 1 

X0UT(I+l)=X0UT(I)+(1.25**(I-l))*RESBL*xvect 
Y0UT(I+1)=Y0UT(I)+(1 . 25**(I-1) ) *RESBL*yvect 
write (* , *) i+1 , X0UT(I+1) ,Y0UT(I+l) 

CONTINUE 

DLEFT= ( 1 . 25** ( 1-1 ) ) * RESBL 



Size the cells at the wall: 

X0UT(NP0INTS)=XEND 
YOUT (NPOINTS ) =YEND 

NOTE: This is a REVERSE counter! ! ! ! ! 

DO 200 1=1 , NNEARWALL- 1 

Each cell is 25*/. larger than the previous cell. 

X0UT(NP0INTS-I )=X0UT (NPOINTS-I+1 )-( 1 . 25** ( 1-1 ) ) *RESWALL*xvect 
Y0UT(NP0INTS-I )=Y0UT (NPOINTS-I+1 ) -( 1 . 25* * ( 1-1 ) ) *RESWALL*yvect 
write ( * , * ) npo int s-i , XOUT (NPOINTS-I ) , YOUT (NPOINTS-I ) 

CONTINUE 

DRIGHT= ( 1 . 25* * ( I- 1 ) ) *RESWALL 
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C Send the remainder of the line to FNSPLT : 

c 

c 

C First need to generate a small line 

c 

c Endpoints of line: XOUT(NNEARBL) , YOUT (NNEARBL) 

C XOUT (NPOINTS-NNEARWALL) , YOUT(NPOINTS-NNEARWALL) 

C FNSPLT WILL NEED 5 POINTS 

TEMPX=XOUT (NPOINTS-NNEARWALL+1 ) -XOUT (NNEARBL ) 
TEMPY=Y0UT(NP0INTS-NNEARWALL+1)-Y0UT(NNEARBL) 
c 

c debugger 

c write (*,*)tempx,tempy, XOUT (NPOINTS-NNEARWALL) 

c 

XLIN( 1 ) =X0UT (NNEARBL) 

YLIN( 1 ) =Y0UT( NNEARBL) 

DO 300 1=1,4 

XLIN( 1+1 ) =XLIN ( I ) +0 . 25+TEMPX 
YLIN( 1+1 ) =YLIN( I ) +0 . 25*TEMPY 
c write(* ,*) ’Straight Line’ 

c write (* , *)xlin(i) ,ylin(i) 

300 CONTINUE 

C 

c 

CALL FNSPLT (5 , NREMAIN , DLEFT , DRIGHT , XLIN , YLIN , XREM , YREM ) 

C 

C FILL OUT THE ARRAY 

DO 400 1=1, NREMAIN 

XOUT ( NNEARBL+I -1 ) =XREM ( I ) 

Y0UT(NNEARBL+I-1 ) =YREM ( I ) 

400 CONTINUE 
C 

C debugger 

c WRITE(* , *) NNEARBL , NNEARWALL , NPO I NTS , NREMA IN 

C 

C debugger 

c 

C 

c write(79,*) ’ZONE’ 

c do 500 i=l,npoints 

c write(79,*) xout (i) , yout (i) 

c 500 continue 
c 

c end debugger 

c 



RETURN 

END 

C+++++++++++++++H^+++^++++++++-H-++++++++++++++++++++++-H-++++++++-H-+++++++ 
SUBROUTINE TUNC (XBEG , XEND , YWALL , NPOINTS , XIN , YIN , XOUT , YOUT ) 

C 

C LAST MODIFIED: 08 NOVEMBER 1996 
C 

C This subroutine takes the cell spacing along the upstream line 

c in the center of the tunnel, the foil points and the down 

c stream line and spaces them along tunnel walls, 
c 

C XBEG = BEGINNING XCOORD OF LINE 
C 

C XEND = ENDING XCOORD OF LINE 
C 

C YWALL = Y VALUE AT WALL 



137 



c 

C NPOINTS = NUMBER OF ARRAY VALUES COMING IN 
C 

C XIN.YIN, = VALUES COMING IN THAT WILL BE PROJECTED 
C 

C XOUT , YOUT= X AND Y VALUES FOR THE LINE ON THE WALL 
C 

IMPLICIT NONE 

REAL P I , TWOP I , ZERO , ONE , HALF 
INTEGER NI,I 

PARAMETER (PI=3. 14159, TW0PI=6. 28318, NI=200, ZER0=0.0, 0NE=1.0) 
INTEGER NPOINTS 

REAL XEND, XBEG , DELX , DELY , TEMP 3 , SLENGTH , TEMP4 , YWALL 
REAL XIN(NI) , YIN(NI) , XOUT(NI), YOUT(NI) .DELARC(NI) 

C 

C Compute arc length of the input array and arc length between 
c points, 
c 

SLENGTH=ZERO 

DELX=ZERO 

DELY=ZERO 

TEMP3=ZER0 

DO 100 1=1, NPOINTS- 1 
DELX=XIN(I+1)-XIN(I) 

DELY=YIN ( 1+ 1 ) -YIN ( I ) 

TEMP3=DELX*DELX+DELY*DELY 
DEL ARC ( I ) =SQRT (TEMP3 ) 

SLENGTH=SLENGTH+DELARC ( I ) 

100 CONTINUE 
C 

C debugger 

C temp4=xin(npoints)-xin( 1) 

C write(* , *)s length , temp4 , ywall , XBEG , XEND 

c 

c do algebraic translation of point spacing 

c 

TEMP4=ZER0 
XOUT ( 1 )=XBEG 
YOUT ( 1 ) =YWALL 
DO 200 I=1,NP0INTS-1 

TEMP4=DELARC ( I ) /SLENGTH 

XOUT ( 1+ 1 )=X0UT ( I ) +TEMP4* ( XEND-XBEG) 

YOUT(I+l)=YWALL 
200 CONTINUE 
C 

C debugger 
C 

c write(79 , *) ’ZONE * 

c do 300 i=l,npoints 

c write(79,*) xout (i) , yout (i) 

c 300 continue 
c 

c end debugger 

c 

RETURN 

END 

SUBROUTINE WALLFIND (SCALE , YUP , YBOT ) 

C 

C LAST MODIFIED: 29 OCTOBER 96 

C 

C This subroutine finds the upper and lower y coordinates of the 

c tunnel walls 
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IMPLICIT NONE 

REAL PI, TWOPI, ZERO, ONE, HALF, TWO 
INTEGER NI 

PARAMETER (PI=3 . 14159 ,TW0PI=6 . 28318 ,NI=10 , 
+ZER0=0 . 0 , ONE= 1.0, TW0=2 . 0 ) 

REAL SCALE, YUP,YB0T 

YUP=0NE/ (SCALE*TW0) 

YB0T=-0NE/(SCALE*TW0) 



debugger 

write(*,*) yup,ybot 



RETURN 

END 

-++++++++++-H-+++-H-++++++++-H-+++-H-+++-H-++-H-+++-H-++++++++-H-+++-H-+++++++ 

SUBROUTINE TUNB (XL , XLOC , YLOC , PLINEX , PLINEY , NP ) 

LAST MODIFIED: 24 OCTOBER 96 

This subroutine defines the upstream or downstream lines extending 
from the foil leading or trailing edge. 

XL = LENGTH OF LINE 

XLOC, YLOC = BEGINNING OF LINE 

PLINEX, PLINEY = OUTPUT X AND Y 

NP = NUMBER OF POINTS OUT 
IMPLICIT NONE 

REAL PI, TWOPI, ZERO, ONE, HALF 
INTEGER NP 

PARAMETER (PI=3 . 14159 , TW0PI=6 . 28318 , ZER0=0.0, 0NE=1.0) 

REAL XL, XLOC, YLOC, TEMP 
INTEGER I 

REAL PLINEX (NP), PLINEY (NP) 



debugger 

write(*,*) ’Made it in to TUNB, LIMIT IS:’, XL 

MAKE IT HORIZONTAL 
DO 100 1=1, NP 
PLINEY ( I ) =YL0C 
100 CONTINUE 

debugger 

write(*,*) xloc,yloc 

The endpoints of the line are xloc,yloc and xloc+temp,yloc 
IF (XLOC. LE. ZERO) THEN 
PLINEX (NP)=XL0C-XL 
PLINEX( 1)=XL0C 
ENDIF 

IF (XLOC. GT. ZERO) THEN 
PLINEX (NP)=XL0C+XL 
PLINEX ( 1 ) =XL0C 
ENDIF 

Generate some points between the two endpoints. 
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TEMP= (PLINEX (NP ) -PLINEX ( 1 ) ) /FLOAT (NP-1 ) 
DO 200 1=1 ,NP-2 

PLINEX (I+1)=PLINEX( I )+TEMP 
200 CONTINUE 

C 

C DEBUGGER 
c write(79,*) 'ZONE’ 

c do 5020 i=l,NP 

c write(79,*) PLINEX(i) ,PLINEY(i) 

c 5020 continue 
RETURN 
END 



C+++++++++++++++++ ++++++++-+-M- ++++++++++++ 

SUBROUTINE ARC JO INER (N 1 , N2 , NOUT , X 1 , X2 , Y 1 , Y2 , XOUT , YOUT ) 

C 

C last modified: 25 October 96 

c 

c This subroutine takes in two arcs and joins them end to end 
c in to one array. 

IMPLICIT NONE 

REAL P I , TWOP I , ZERO , ONE , HALF 
INTEGER NI , I 

PARAMETER (PI=3 . 14159 ,TW0PI=6 . 28318 , NI=200, ZER0=0.0, 0NE=1.0) 
INTEGER Nl, N2,N0UT 

REAL XI (NI) , Yl(NI), XOUT(NI), YOUT(NI), X2(NI) , Y2(NI) 

C 

C debugger 

c write(*,*) 'Made it in to ARCJOINER: > , nl,n2 

c 

c 

DO 10 1=1, Nl 
X0UT(I)=X1(I) 

Y0UT(I)=Y1(I) 

10 CONTINUE 

DO 20 1=2, N2 

XOUT ( I+Nl-1 )=X2 ( I ) 

Y0UT(I+N1-1)=Y2(I) 

20 CONTINUE 

c 

c funny counting here elmininates overlapping data points 

c 

N0UT=N1+N2-1 

C 

C debugger 

c write(*,*) 'Points Out',nout 

c 

c 

RETURN 

END 

C+++++++++++++++-H-+++-H-++++++++-H-+++-H-+++-H-+++++++++++++++++++++++++++++++ 

SUBROUTINE ARCMAKER ( X , Y , NUM 1 , NUM2 , XO , YO ) 

C 

C last modified: 29Sept96 

C 

C THIS SUBROUTINE TAKES IN LONG ARRAYS X, AND Y, 

C AND RETURNS PORTIONS OF X AND Y SPECIFIED BY 

C INTEGERS NUM1 AND NUM2 . XO AND YO ARE FILLED 

C BY THE FIRST NUM1-NUM2 ELEMENTS 

C 

IMPLICIT NONE 
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REAL P I , TWQPI , ZERO , ONE , HALF 
INTEGER NI,I 

PARAMETER (PI=3. 14159, TW0PI=6. 28318, NI=200, ZER0=0.0, 0NE=1.0) 

INTEGER NUM1.NUM2 

REAL X(NI) , Y(NI),XO(NI), YO(NI) 

C 

C debugger 
c 

C write(* , *)NUM1 ,NUM2 

DO 10 1=1 , (NUM2-NUM1+1 ) 

X0(I)=X(NUM1+I-1) 

Y0 (I)=Y (NUM1+I-1 ) 

debugger 

write(* ,*) X0(I),Y0(I) 

10 CONTINUE 
RETURN 
END 

SUBROUTINE F0ILZ0NE(XI ,NP0INTS ,MIDPRES , MIDSUCT ,N0SE) 
last modified: 29 Sept 96 

THIS SUBROUTINE LOOKS AT ALL THE FOIL GEOMETRY POINTS. IT SEARCHES 
FOR THE LEADING EDGE OF THE FOIL AT THE SET AOA. IT SEARCHES FOR 
POINTS ON THE SUCTION AND PRESSURE SIDES THAT ARE CLOSEST TO THE 
MIDCHORD POINT. ONCE THESE POINTS ARE IDENTIFIED, THE INTEGER 
OF THE ARRAY POSITION OF THESE POINTS IS RETURNED TO MAIN PROG 



MIDSUCT 

I 

NOSE < ================== ==\ (foil schematic) 

I 

MIDPRES 



IMPLICIT NONE 

REAL PI , TWOPI , ZERO , ONE , HALF 
INTEGER NI , I 

PARAMETER (PI=3 . 14159 ,TV0PI=6 . 28318 , NI=200, ZER0=0.0, 0NE=1.0) 
INTEGER NPO INTS , MIDPRES , MIDSUCT , NOSE , TEMP 
REAL CHORDMID 
REAL XI (NI) 

XI, YI ARE FOIL OFFSETS 

NPOINTS : TOTAL NUMBER OF FOIL OFFSETS 

MIDPRES, MIDSUCT AND NOSE AS INDICATED ON SCHEMATIC 

FIND THE LEADING EDGE: 

TEMP=0 

DO 10 1=1 , ( NPOINTS- 1) 

IF ( (XI (I) ) . GE . (XI (I+l) ) ) THEN 
TEMP=I+1 
ENDIF 

10 CONTINUE 
NOSE=TEMP 

FIND THE MIDDLE OF THE PRESSURE SIDE 
TEMP=0 

CHORDMID= ( ( X I ( NOSE) +XI ( 1 ) ) ) / ( 2 . 0 ) 
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chordmid 



c 

c debugger 

c write(*,*) 'pressmid= J , 

c 

DO 20 1=1, (NOSE-1) 

IF (XI ( I ).GT. CHORDMID) THEN 
TEMP=I 
ENDIF 

20 CONTINUE 

MIDPRES=TEMP 

C 

C FIND THE MIDDLE OF THE SUCTION SIDE 

c 

TEMP=0 

CHORDMID=(XI (NOSE)+XI(NPOINTS) )/(2 . 0) 
c write(*,*) ’ suctmid^, chordmid 

DO 30 I=NOSE, (NPOINTS-1) 

IF (XI ( I ).LT. CHORDMID) THEN 
TEMP=I 
ENDIF 

30 CONTINUE 

MIDSUCT=TEMP 

c 

c debugger 
c 

c WRITE ( * , * ) NPOINTS , NOSE , MIDSUCT , MIDPRES 

c 

RETURN 

END 

C 

C+++++++++++++++-H-++++++++-H-++++++++-H-++++++++++++-H-++++++++++++++-f+++++++ 
SUBROUTINE NACACONV (N , X , Y , T , R) 

C 

C last modified: 23 Sept 96 

C 

C THIS SUBROUTINE CONVERTS NACA FOIL GEOMETRY TO X,Y GEOMETRY 
C 

C N = NUMBER OF STATIONS 

C X = STATION 

C Y = CAMBER 

C T = THICKNESS 

C R = L.E. RADIUS 

C 

C RETURNS 

C N = NUMBER OF POINTS (.GT. 2*STATI0NS) 

C X , Y = GEOMETRY COORDINATES 

C 

IMPLICIT NONE 

REAL P I , TWOP I , ZERO , ONE , HALF 
INTEGER NI 

PARAMETER (PI=3 . 14159, TW0PI=6 . 28318 , NI=200) 

INTEGER I ,N 
REAL R 

REAL X(NI),Y(NI),T(NI) 

THIS SUBROUTINE NOT YET FULLY IMPLEMENTED 

WRITER,*) * THE NACA CONVERSION ROUTINE IS NOT YET IMPLEMENTED * 
WRITE(* , * ) ’ Pretty lame, huh?’ 

END 

+++++++++++++++^+++H^+++4H-+++-^+++-H-+++++++-H-+++-H-+++-M-+++++++++++++-H*++ 
SUBROUTINE ROTANGLE (N , X , Y , AOA , PIVOTX , PIVOTY) 

last modified: 26 Sept 96 
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c 

C THIS SUBROUTINE TAKES X,Y FOIL DATAPOINTS AND ROTATES THEM ABOUT 

C A PIVOT POINT TO SET FOIL TO PRESCRIBED ANGLE OF ATTACK(AOA) 

C 

C N = NUMBER OF POINTS DEFINING FOIL 
C X = X COORDINATE OF POINT (X/C) 

C Y = Y COORDINATE OF POINT (Y/C) 

C AOA = FOIN ANGLE OF ATTACK IN DEGREES REF TO NOSE TAIL LINE 

C PIVOT = X/C OF ROTATION POINT OF FOIL 

C 

C XL, YL ARE REF TO PIVOT POINT 
C R IS RADIUS TO PIVOT POINT 
C 

IMPLICIT NONE 

REAL PI, TWOPI, ZERO, ONE, HALF 
INTEGER NI 

PARAMETER (PI=3 . 14159, TW0PI=6 . 28318, NI=200) 

INTEGER I ,N 

REAL AOA, PIVOTX , PIVOTY, XL, YL, R, BETA 
REAL X(NI),Y(NI) 

C 

C 

WRITE(* , *) 1 FOIL BEING ROTATED TO ANGLE OF ATTACK: 1 , AOA 
WRITE(* , *) ' PIVOT POINT IS X COORD <Q (X/C) : 1 , PIVOTX 
WRITE(* , *) 1 PIVOT POINT IS Y COORD <0 (Y/C) : 1 , PIVOTY 

C 

C ROTATE FOIL POINT BY POINT 
C 

DO 101 I =1 ,N 
XL=X(I) -PIVOTX 
YL=Y ( I ) -PIVOTY 
R = SQRT(XL*XL+YL*YL) 

BETA = ATAN2D(YL,XL) 

X(I) = R*COSD(BETA-AOA) 

Y(I) = R*SIND(BETA-AOA) 

101 CONTINUE 
RETURN 
END 
C 

C 

SUBROUTINE NORMFOIL(L,N,X,Y) 
c 

c LAST MODIFIED: 24 OCTOBER 96 

IMPLICIT NONE 

REAL PI, TWOPI, ZERO, ONE, HALF 
INTEGER NI , I 

PARAMETER (PI=3 . 14159 , TW0PI=6 . 28318 , NI=200, ZER0=0.0, 0NE=1.0) 
INTEGER N 
REAL L 

REAL X(NI) , Y(NI) 

DO 10 1=1, N 
X(I)=X(I)/L 
Y(I)=Y(I)/L 
10 CONTINUE 
RETURN 
END 



C+++++++++++++++-4-f+++++++++++++++++++++++++++++++++++HH-++++++++++++++++++ 

SUBROUTINE FNSPLT (NIN , NOUT , DS 1 , DS2 , XI , YI , XO , YO ) 

C 

C | i m | m | | j j | m i j i m i j j i ! i j | | m j j j J m j ! ! ! ! j j j ! J j ! i ! i j J j j i j | j j m j j 

C ! THIS PROGRAM* CONVERTED IN TO A SUBROUTINE IN M FIT2D . F" ! 
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c ! BY JOHN DANNECKER. THIS PROGRAM IS ORIGINALLY WRITTEN ! 

C ! BY S.D. BLACK, MIT, WITH OTHER CREDITS AS INDICATED ! 

C !!!!!!!!!!! H !!!!!!■!!! M MM 1 1 i i i i i i i i i i i i M i i i i i i i i ii i i i i i i i i i 

C 

C FANCY SPLITTING PROGRAM FOR DIVIDING UP A CURVE 
C THERE ARE TWO METHODS THIS PROGRAM CAN BE RUN 
C 

C INPUT 1 : 

C DATA SET (X,Y) FOR THE CURVE 

C DS1 - INTERVAL AT LEFT HAND END POINT 

C DS2 - RIGHT HAND END POINT INTERVAL 

C SET NOUT < 0 

C 

C INPUT 2 : 

C SAME AS INPUT 1 BUT WITH THE NUMBER OF POINTS SPECIFIED 

C 

C 

C OUTPUT 

C X,Y OF SPLIT UP CURVE BETWEEN THESE INTERVALS 

C 

C PURPOSE 

C USEFUL FOR CUSTOM GRIDDING OF CURVES 

C 

C THE SLOPE OF THE INTERVAL DS# IS HELD TO BE ZERO AT THE END POINTS 

C THE PROGRAM ASSUMES THE ENDS OF THE CURVE ARE AT THE ENDS OF THE INPUT 

C 

C WRITTEN 1/8/95 S. BLACK, MOD BY J. DANNECKER 9/28/96 
C 

c This is the NEW and Improved FNSPLT modified by 
c S. Black 25 October 1996 
c 

PARAMETER (NI=200) 

REAL XI (NI) ,YI(NI) , XO(NI) , YO(NI) , CUB(4*(NI-1) ) 

REAL SI(NI) , CUB1(4* (NI-1) ) , CUB2(4* (NI-1) ) 

REAL A(3),B(3),C(NI),D(NI),E(NI) 

CHARACTER FIN*20 
CHARACTER PR0MPT2+34 



debugger 

write(* ,*) ’made it in to f nsplt * , 
################################## 



nin,nout 



CALCULATE LENGTH OF CURVE BASED ON INPUT POINTS 
SL=0 . C 

DO 10 1=2, NIN 

SL=SL+SQRT ( (XI (I)-XI (1-1) ) **2+(YI (I) -YI (1-1) ) **2) 

10 CONTINUE 

CALCULATE THE NUMBER OF POINTS REQ'D IF NOT SPECIFIED 
IF (NOUT .LE . 0) THEN 
DSAVG=(DSl+DS2)/2 . 0 
N0UT=INT(SL/DSAVG)+1 
END IF 

SPLINE INPUT ARRAY PARAMETRICALLY (BOTH X AND Y) 

Array SI is returned with arclength parameters non-dimed from 0..1 
CALL PUGLYDK (NIN , XI , YI , SI , CUB 1 , CUB2 ) 

SET UP ARRAYS CONTAINING SPACING 

Array A Contains pointers for the first, middle and last 
steps . 

Array B contains the dS values at the two ends and a guess 
at what the step size in the middle should be. 

Array C contains integers cooresponding to the NOUT-1 steps 
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DSAVG= ( SL-DS 1-DS2 ) /FLO AT (NOUT-3 ) 

A ( 1 ) =FLOAT ( 1 ) 

A ( 2 ) =FLOAT (NOUT ) /2 . 0 
A(3)=FL0AT( NOUT- 1 ) 

B(1)=DS1 

B(2)=DSAVG 

B(3)=DS2 

DO 20 J=1,N0UT-1 
C(J)=FLOAT(J) 

20 CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

DSMIN=0.1*MIN(DS1,DS2) 

C. . .CONVERGENCE LOOP FOR APPROPRIATE VALUE OF B(2) 

C The step size is splined with the slope held constant at the 
C ends. There must be NOUT-1 steps, with values of DS1 and 
C DS2 at the two ends. The mid-range step size is varied until 

C the sum of all the steps equals the arc length desired (SL) . 

C By fixing the slope at the ends the program tries to ensure 

C that step sizes don’t increase too rapidly. 

DO 100 KK=1 , 1000 
IERR=0 

C... SPLINE SPACING WITH FIXED SLOPE AT THE ENDS 
CALL UGLYDK (3,0,0,A,B,0,0, CUB ) 

C... EVALUATE SPLINE TO FIND STEP SIZE AT INTERMEDIATE POINTS 
CALL EVALDK (3 , NOUT-1 ,A,C,D,CUB) 

C. . .CALCULATE LENGTH COVERED BY NEW SET OF STEPS 
SLC=0.0 

DO 40 1=1, NOUT-1 
SLC=SLC+D(I) 

IF (D(I) .LT.O.O) IERR=1 
40 CONTINUE 

C.. .IF LENGTH IS NOT CLOSE TO TOTAL LENGTH NEEDED, SHIFT B(2) 

C AND REITERATE 

IF (IERR.EQ.l) THEN 
B(2)=0.0 

ELSE IF (ABS(SLC-SL) .GT.DSMIN) THEN 
DIFF=SLC-SL 
B(2)=B(2)-DIFF/N0UT 
ELSE 

GOTO 101 
END IF 

100 CONTINUE 

WRITE(* ,*) ’ ***WARNING*** FNSPLT DID NOT CONVERGE’ 

WRITE(*,*) ’ ABS(SLC - SL) > TOL’ 

WRITE(*,*) ’ SLC , SL, TOL = ’ ,SLC,SL,DSMIN 
WRITE(*,*)NIN,N0UT,DS1,DS2 
IF (SLC.GT.SL) THEN 

WRITE(*,*) ’Recommend using fewer points or smaller dS values’ 
ELSE 

WRITE(*,*) ’Recommend using more points of larger dS values’ 
END IF 

101 CONTINUE 

ccccccccccccccccccccccccccccccccccccccccccccccccc 
C... SHIFT LENGTH TO EXTEND FROM 0.0 TO 1.0 
E(l)=0.0 
DO 50 1=2, NOUT-1 

E(I)=D(I-1)/SLC+E(I-1) 

50 CONTINUE 

E(N0UT)=1 . 0 

C... EVALUATE OUTPUT POINTS AT CALCULATED INTERVAL 
CALL PEVALDK (NIN , NOUT , E , SI , XO , YO , CUB 1 , CUB2 ) 

C 

c 

C ################################### 
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cn cr> 



RETURN 

END 



SUBROUTINE DRI VDK (NIN , NOUT , X IN , XOUT , D YDX , D2YDX , A ) 

C APRIL 1975 SPLINE PROGRAM SERIES J.E.KERWIN 
REAL XIN(*) ,XOUT(*) , DYDX(*) ,D2YDX(* ) , A(*) 

NM1=NIN-1 

J=1 

DO 3 N=l,NOUT 

IF(X0UT(N).GE.XIN(2)) GO TO 4 
J=1 

GO TO 5 

4 IF(XOUT(N) .LT.XIN(NMl) ) GO TO 6 
J=NM1 
GO TO 5 

IF(XOUT(N).GE.XIN(J+l)) GO TO 7 
Hl=XOUT(N)-XIN(J) 

H2=H1**2 
J2=J+NM1 
J3=J2+NM1 

DYDX(N)=3 . 0*A( J)*H2+2 . 0*A(J2) *H1+A( J3) 

D2YDX(N)=6 . 0*A( J) *Hl+2 . 0*A( J2) 

GO TO 3 
7 J=J+l 

GO TO 6 

3 CONTINUE 
RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE EVALDK (NIN , NOUT , XIN , XOUT , YOUT , A ) 

C APRIL 1975 SPLINE PROGRAM SERIES J.E.KERWIN 
REAL XIN(*) ,XOUT(*) ,YOUT(*) , A(*) 

NM1=NIN-1 
MOUT=I ABS (NOUT ) 

IF (NOUT . GT . 0) GO TO 1 
DEL=(XIN(NIN)-XIN( 1 ) )/ (MOUT-1 ) 

DO 2 N=1 ,MOUT 

2 X0UT(N)=XIN(1)+(N-1)*DEL 

1 J=1 

DO 3 N=1,M0UT 

IF(XOUT(N) .GE.XIN(2)) GO TO 4 
J=1 

GO TO 5 

4 IF ( XOUT (N) .LT.XIN(NMl) ) GO TO 6 
J=NM1 

GO TO 5 

6 IF(XOUT(N) .GE.XIN(J+l) ) GO TO 7 
9 IF (XOUT(N) .LT.XIN(J) ) GO TO 8 

5 H1=X0UT(N)-XIN(J) 

H2=Hl**2 

H3=H1*H2 

J2=J+NM1 

J3=J2+NM1 

J4=J3+NM1 

YOUT (N)=A( J) *H3+A( J2) *H2+A( J3) *H1+A( J4) 

GO TO 3 

7 J=J+1 
GO TO 6 

8 J=J-1 
GO TO 9 

3 CONTINUE 
RETURN 
END 

C 

C+++++++++++++++-H-+++-H-+++-H-++++++++-H-++++++++++++++++++++++++++++++++++++ 
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C 

SUBROUTINE INTEDK (NIN , XIN , XL , XU , YDX , XYDX , XXYDX , A ) 

C AUGUST 1975 SPLINE PROGRAM SERIES S.-K.TSAO 

REAL XIN(*),A(*) 

NM1=NIN-1 

IF(XL.LT. XIN(D) GO TO 2 
DO 1 N=1 ,NIN 
IF(XL.GE . XIN(N) ) GO TO 1 
JL=N-1 
GO TO 3 
1 CONTINUE 

J L=NM 1 
GO TO 3 
JL=1 

IF (XU. GE. XIN (NIN) ) GO TO 4 
JU=1 

IF (XU . LE . XIN ( 1 ) ) GO TO 6 
DO 5 N=JL,NIN 
IF(XU .GE. XIN(N) ) GO TO 5 
JU=N-1 
GO TO 6 

5 CONTINUE 
GO TO 6 

4 JU=NM1 

6 H1=XL-XIN(JL) 

H2=H1**2 

H3=H1*H2 

H4=H2**2 

H5=H2*H3 

H6=H3**2 

J 1=JL 
J2=J 1+NM1 
J3=J2+NM1 
J4=J3+NM1 

YDX=-A( Jl)/4 . 0*H4-A( J2)/3 . 0*H3-A( J3)/2 . 0*H2-A( J4) *H1 
BUG=-A( Jl)/5 . 0*H5-A( J2)/4 . 0*H4-A( J3)/3 . 0*H3-A( J4)/2 . 0*H2 
XYDX=BUG+XIN ( J 1 ) * YDX 

BUG=-A( Jl)/6 . 0*H6-A( J2)/5 . 0*H5-A ( J3)/4 . 0*H4-A( J4)/3 . 0*H3 
XXYDX=BUG+2 . 0*XIN( Jl) +XYDX-XIN ( J 1 ) * *2*YDX 
DO 7 N=JL , JU 
H1=XIN(N+1)-XIN(N) 

IF(N.EQ.JU) H1=XU-XIN(N) 

H2=H1**2 
H3=H1*H2 
H4=H2**2 
H5=H2*H3 
H6=H3**2 
J 1=N 

J2=J 1+NM1 
J3=J2+NM1 
J4=J3+NM1 

BUG=A(J1)/4.0*H4+A(J2)/3.0*H3+A(J3)/2.0*H2+A(J4)*H1 

YDX=YDX+BUG 

CAT=A(Jl)/5 .0*H5+A( J2)/4 .0*H4+A( J3)/3 .0*H3+A(J4)/2. 0*H2 

PIG=CAT+XIN( J l) *BUG 

XYDX=XYDX+PIG 

D0G=A(J1)/6.0*H6+A(J2)/5.0*H5+A(J3)/4.0*H4+A(J4)/3.0*H3 
XXYDX=XXYDX+D0G+2 . 0*XIN( Jl) *PIG-XIN( Jl) **2*BUG 

7 CONTINUE 
RETURN 
END 

C 

0 +++++++++++++++^++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

C 

SUBROUTINE LINDK(NIN,XIN, YIN, A) 



147 



o o 



C GENERATES CUBIC COEFFICIENTS FOR PIECEWISE LINEAR FIT 

REAL XIN(*),YIN(*),A(*) 

NM1=NIN-1 
DO 1 N=1 ,NM1 
A(N)=0 . 0 
M=N+NM1 
A(H)=0 . 0 
M=M+NM 1 

A(M)=(YIN(N+1)-YIN(N))/(XIN(N+1)-XIN(N)) 

M=M+NM1 

1 A(M)=YIN(N) 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE PEVALDK (NIN , NOUT , ARK , S , XO , YO , CUB 1 , CUB2) 

REAL X0(200) ,Y0(200) ,S(200) ,ARK(200) 

REAL CUB1( (200-1) *4) ,CUB2( (200-1) *4) 

CALL EVALDK (NIN , NOUT , S , ARK , XO , CUB 1 ) 

CALL EVALDK (NIN , NOUT , S , ARK , YO , CUB2) 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE PUGLYDK(NIN,X,Y,S,CUB1,CUB2) 

REAL X(200) ,Y(200) ,S(200) 

REAL CUB1((200-1)*4),CUB2((200-1)*4) 

ST0T=0 . 0 
S(1)=0.0 
DO 10 1=2, NIN 

S(I)=SQRT( (X(I)-X(I-l) )**2+(Y(I)-Y(I-l) )**2)+ST0T 
STOT=STOT+SQRT( (X(I)-X(I-l) )**2+(Y(I)-Y(I-l) ) **2) 

10 CONTINUE 

DO 20 1=2, NIN 

S(I)=S(I)/STOT 
20 CONTINUE 
NCL=1 
NCR=1 
ESL=0 
ESR=0 

CALL UGLYDK (N IN , NCL , NCR , S , X , ESL , ESR , CUB 1 ) 

CALL UGLYDK (NIN , NCL , NCR , S , Y , ESL , ESR , CUB 2) 

RETURN 

END 

C 

SUBROUTINE UGLYDK (NIN , NCL , NCR , XIN , YIN , ESL , ESR , AE ) 

-1975 DUCK SERIES J.E.KERWIN MODIFIED 6/21/82 

-TRI-DIAGONAL MATRIX SOULTION BUILT IN 

REAL XIN(*) ,YIN(*) , AE(*) ,H(200) ,D(200) ,AU(200) ,AM(200) , 

* S(200) ,AL(200) ,X(200) 

DATA HALF/O . SEOO/ , TWO/2 . OEOO/ , SIX/6 . OEOO/ , RAD/1 . 745329E-02/ 
NM1=NIN-1 
NM2=NM1-1 
NM3=NM2-1 
NEQ=NM2 
DO 1 N=1 ,NM1 
H(N)=XIN(N+1)-XIN(N) 

1 D(N)=(YIN(N+1)-YIN(N))/H(N) 

IF(NCL.EQ.2) NEq=NEQ+l 
IF (NCR. EQ. 2) NEq=NEQ+l 
NSq=NEq**2 
J=1 

IF(NCL. LT. 2) GO TO 6 
AM(l)=TWO*H(l) 

AU ( 1 ) =H ( 1 ) 

SLP=ESL*RAD 

S(1)=(D( l)-TAN(SLP) ) *SIX 
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J=J+1 

AL(2)=H( 1) 

DO 5 N=1,NM2 
IF(N.GT.l) AU( J-1)=H(N) 

AM(J)=TWO*(H(N)+H(N+l) ) 

IF(N.LT.NM2) AL( J+1)=H(N+1) 

IF(N.EQ . 2. AND.NCL. Eq . 1) AU(J-1)=AU(J-1)-H(N-1)**2/H(N) 

IF (N . EQ . 1 . AND . NCL . EQ . 1 ) AM( J)=AM( J)+( 1 . 0+H(N)/H(N+l) ) *H(N) 
IF(N.Eq.NM2.AND.NCR.EQ.l) AM(J)=AM(J)+(1.0+H(N+l)/H(N))*H(N+l) 
IF(N.EQ.NM3.AND.NCR.EQ. 1) AL ( J+l ) =AL ( J+ 1 ) -H (N+2 ) * *2/H(N+l ) 
S(J)=(D(N+1)-D(N) )*SIX 
J=J+1 
5 CONTINUE 

IF(NCR.LT. 2) GO TO 7 
AL(NEQ)=-H(NM1) 

AM(NEq)=-TWO*H(NMl) 

AU(NEq-l)=H(NMl) 

SLP=ESR*RAD 

S( J)=(D(NM1)+TAN(SLP) ) *SIX 

7 CONTINUE 

DO 4 K=2,NEq 
AL(K)=AL(K)/AM(K-1) 

AM(K)=AM(K)-AL(K) *AU(K-1) 

S(K)=S(K)-AL(K)*S(K-1) 

4 CONTINUE 

X(NEq) =S(NEq ) /AM (NEq ) 

DO 2 L=2,NEq 
K=NEq-L+l 

X(K)=(S(K)-AU(K)*X(K+1))/AM(K) 

2 CONTINUE 

DO 22 N= 1 ,NEq 
22 S(N)=X(N) 

HOLD=S(NEq) 

IF(NCL.Eq.2) GO TO 8 
DO 9 N=1,NM2 
M=NM2-N+2 

9 S(M)=S(M-1) 

IF(NCL.Eq.O) S( 1 )= 0.0 
BUG=H( 1)/H(2) 

IF(NCL.Eq.l) S(l)=( 1 . O+BUG) *S(2)-BUG*S(3) 

8 IF(NCR.Eq.O) S(NIN)=0 . 0 
BUG=H (NM1)/H(NM2) 

IF(NCR.Eq.l) S(NIN)=( 1 . O+BUG )*S(NM1)-BUG*S(NM2) 

IF(NCR.Eq.2) S(NIN)=HOLD 
DO 10 N=1 ,NM1 

AE(N)=(S(N+1)-S(N) )/ (SIX*H(N) ) 

M=N+NM 1 

AE(M)=HALF*S(N) 

M=M+NM1 

AE(M)=D(N)-H(N) *(TU0*S(N)+S(N+1) )/SIX 
M=M+NM 1 

10 AE(M)=YIN(N) 

RETURN 

END 
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Appendix B 

Sample FIT2D Input Files 



B.l Bounded Foil (fname. Ctrl) 

Header: Sample for Cupped B-l 

Foil geometry file name 

rcup.foil 

AOA Xpiv Ypiv 

0.5 0.3 0.01 

Chord/Tunnel Width 

0.9 

USL DSL PHI1 PHI2 

2.0 2.0 0.5 0.25 

NUS NDS NVS NT0P NB0T (Always 8*n-l points) 
87 87 95 95 95 

RESLE RESMID RESTE RESWAL RESBL PACK 
0.001 0.03 0.001 1.0e-5 1.0e-5 0.20 

Re#(chord based) 

3e6 

Desired Inmesh Input file: 
cupi . dat 

Desired Inmesh Restart file 
cupr . dat 

Number of INMESH Iterations: 

1500 

Convergence tolerance: 

1.0e-10 

Wake Data(-1 none, 0 VLM, or nranswk) 

10 

0.713324 -0.016042 
0.772837 -0.022880 
0.879395 -0.032600 
1.027451 -0.043110 
1.210484 -0.052885 
1.424693 -0.061495 
1.674272 -0.069320 
1.953325 -0.077086 
2.257376 -0.083998 
2.580143 -0.089997 
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B.2 Unbounded Foil (/name. Ctrl) 

Header: Sample Unbounded Cupped B-l Foil 
Foil geometry filename: 
rcup .f oil 
AOA Xpiv Ypiv 
0.5 0.3 0.01 

Chord/Tunnel Width 
0.08 

USL DSL PHI1 PHI2 

5.0 5.0 2.00 1.00 

NUS NDS NVS NT0P NB0T(Always 8*n-l points) 
87 87 95 95 95 

RESLE RESMID RESTE RESWAL RESBL PACK 
0.001 0.03 0.001 0.005 1.0e-5 0.20 

Re# (chord based) 

3e6 

Desired Inmesh Input file: 
cupi.dat 

Desired Inmesh Restart file 
cupr . dat 

Number of INMESH Iterations: 

1500 

Convergence tolerance: 

1.0e-10 

Wake Data(-1 none, 0 VLM, or nranswk) 

20 

0.713324 -0.016042 
0.772837 -0.022880 
0.879395 -0.032600 
1.027451 -0.043110 
1.210484 -0.052885 
1.424693 -0.061495 
1.674272 -0.069320 
1.953325 -0.077086 
2.257376 -0.083998 
2.580143 -0.089997 
2.917188 -0.095028 
3.261140 -0.099227 
3.603834 -0.102689 
3.939133 -0.105482 
4.261238 -0.107661 
4.565246 -0.109321 
4.846685 -0.110522 
5.101448 -0.111331 
5.326634 -0.111836 
5.518953 -0.112118 
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B.3 Sample Foil Geometry File (/name. foil) 



Offsets start at the trailing edge marching forward on lowersurface, around the leading 
edge and along upper surface to trailing edge. It is not neccessary to have an explicit 
point at the leading edge. FIT2D splines the offsets and finds the leading edge. 



FIT2D assumes that this input file is at zero degrees angle of attack. 

Cupped B-l Foil 
2 

1.0 

161 

1.000000 0.000000 

0.989517 0.002052 

0.979110 0.003558 

0.968851 0.004586 

0.958810 0.005204 

0.949059 0.005478 

0.939668 0.005477 

0.930708 0.005269 

0.922250 0.004920 

0.914365 0.004499 



0.007665 - 0.006030 
0.005500 - 0.005229 
0.003652 - 0.004320 
0.002127 - 0.003289 
0.000963 - 0.002132 
0.000220 - 0.000846 
0.000041 0.000570 
0.000133 0.001781 
0.000676 0.003074 
0.001617 0.004450 
0.002986 0.005909 
0.004814 0.007453 
0.007129 0.009082 
0.009962 0.010796 



0.953672 0.028109 
0.962481 0.024014 
0.971452 0.019244 
0.980651 0.013711 
0.990145 0.007326 
1.000000 0.000000 
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Appendix C 

Marine Hydrodynamics Lab Water 
Tunnel Geometry 

C.l System Overview 
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Appendix D 



GETWAKE Program Listing 



PROGRAM GETWAKE 
c 
c 

c Reads in a set of wakepoints . 

INPUT REQUIREMENTS: 

JUNK header line 
NPOINTS - NUMBER OF OFFSETS 
XIN, YIN 
• > 

(to end of file) 
c 
c 
c 

IMPLICIT NONE 

REAL XIN(IOOO) ,YIN(1000) 

REAL Xout(200) ,Yout(200) 

integer npoints ,nf oil , i ,k,ntemp,nout , j , ierr 

character FN0PEN*20 

character junk*30, title*30 , JUNK2*2 

character PR0MPT2*30 

character FIN*20 

WRITE(*,*) ’***’ 

write(*, ’ (A,$) ’ ) ’Enter Filename of Tecplot stream data: ’ 
read(* , ’ (A) ’ ) FNOPEN 
WRITE(* , *) ’***’ 

write(* , ’ (A) ’ ) ’Shortening file:’ // FNOPEN 

OPEN (UNIT=1 , FILE=FNOPEN , F0RM= ’ FORMATTED ’ , STATUS= ’ OLD ’ ) 
WRITE(* , *) ’ *** ’ 

WRITE(* , *) ’ ***READING DATA IN*** ’ 

WRITE(* , *) ’ *** ’ 
read ( 1 , ’ (A) ’ ) junk 
read(l , ’ (A) ’ ) junk 
read(l, ’ (A) ’ ) junk 

C WRITE (*,*) ’READ THE TOP OF THE FILE’ 

c 

c keep next read as the title 
read( 1 , ’ (A) ’ ) title 
c 

c now want to extract the integer for number 
c of pairs of points 

CALL READIJ(1, npoints, J,IERR) 

98 FORMAT (A23 , 16 , A5) 

WRITE(* ,98) ’ ***NUMBER DATAPOINTS IN: ’, NPOINTS,’ ***’ 
WRITE(* , *) ’***’ 
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READ(1, '(A) ') JUNK 

READ(1,*) ( (xin(k) ,yin(k) ) ,K=1 ,npoints) 

C READ(1,*) JUNK 

CLOSE (1) 

WRITE(* , * (A ) * ) 9 ***CLOSING FILE: 1 // FNOPEN // ' ***' 
WRITE(* ,*)’*** » 

c 

c Shorten Data file by taking only 20 datapoints in range 
c 

N0UT=20 

WRITE(*,98) '***NUMBER OF POINTS OUT: ' ,N0UT, 1 ***> 
WRITE (* , * ) ' *** ' 

NTEMP=int (float (npo int s ) /NOUT) 

WRITE(*,98) >***DATA SKIP INTERVAL: ' ,NTEMP, ' ***' 
WRITE(* , *) ' *** 1 

do 100 i=l ,N0UT 

xout ( i) =xin ( i*ntemp-ntemp+l ) 
y ou t ( i ) =y in ( i *nt emp-nt emp+ 1 ) 

100 continue 

write(* ,*) '****Data has been shortened.****** 1 
write(* ,*) '*** ' 



PR0MPT2 = 'Output FILE V/ * ('// 'wakelin.dat '//' ) = 1 
WRITE (* , ' (A , $) ' ) PR0MPT2 
READ (*,'(A)') FIN 

IF (FIN(1:1).EQ. ' ') FIN = 'wakelin.dat' 

WRITE (* , ' (A) ' ) 'OPENING FILE: ' // FIN 

OPEN (UNIT=3 , FILE=FIN , STATUS= ' UNKNOWN ' ) 

91 format (i4) 

92 format (2f 10.6) 

WRITE (*, *) '*** ' 

WRITE(* , *) ' *** WRITING OUTPUT TO FILE***' 

WRITE(* , *) '***' 

WRITE (3,*) 'Shortened: ', TITLE 
WRITE (3 , 91)N0UT 
do 200 i=l,nout 

write(3,92) xout (i) ,yout (i) 

200 continue 
close(3) 

write(*,*)' All Done, Thanks. ' 

WRITE(*,*) '***' 
stop 
end 

cXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

c 

c 

*************************** SUBROUTINE READIJ ********************* 

* EXTRACTS THE INTEGER I,J INDICES FROM A TECPLOT FILE HEADER * 

* Arguments: IUNIT . . . .Unit number of Tecplot file * 

* I Number of columns in IJ ordered data * 

* J Number of rows in IJ ordered data * 

* I ERR 0=Valid data, l=Subroutine failed * 

* It is assumed that I and J are between 1 and 9999 * 

* Justin E. Kerwin August 19,1993 * 

******************************************************************* 

SUBROUTINE READIJ (IUNIT, I , J, IERR) 

CHARACTER*80 LABEL 
CHARACTER*4 IC0DE 
IC0DE=' * 
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IERR=0 

C Read in the Tecplot header line as a character string LABEL. 

READ ( IUNIT , * ( A ) * ) LABEL 

C Find the length of the character string 

LMAX=LEN(LABEL ) 

LBEGIN=1 

C First find I:(IJ=1), then find J:(IJ=2) 

DO 400 IJ=1 ,2 
IC0DE => 

C Find the position of the first = sign in the string 

DO 100 L=LBEGIN,LMAX 

IF(LABEL(L:L).Eq.CHAR( 61 )) THEN 
LMIN=L +1 
GO TO 110 
END IF 



100 CONTINUE 

IERR=1 
GO TO 9999 
110 CONTINUE 

C Find the position of the first number following am = sign. 



DO 200 L=LMIN,LMAX 

IF(LABEL(L:L) .GT. CHAR (48) .AND .LABEL (L :L) . LT. CHAR(58) ) THEN 
LSTART=L 
GO TO 210 
END IF 



200 CONTINUE 

IERR=1 
GO TO 9999 
210 K=1 

C Generate a substring ICODE consisting of consecutive numbers 



DO 300 L=LSTART , LSTART+4 

IF(LABEL(L:L).LT. CHAR ( 48 ) . OR . LAB EL ( L : L ) . GT . CHAR ( 57 ) ) THEN 
LBEGIN=L 
GO TO 310 
ELSE 

ICODE(K:K)=LABEL(L:L) 

K=K+1 
END IF 



300 CONTINUE 

IERR=1 
GO TO 9999 

C Convert the substring to integers I,J using an internal read 

310 IF(IJ.EQ.I) READ (ICODE, * (14) * ) I 

IF(IJ .EQ .2) READ ( ICODE, >( 14)0 J 



400 CONTINUE 

9999 RETURN 
END 
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Appendix E 



PATCH Progam Listing 



PROGRAM PATCH 

REAL*8 X(200 , 200) , Y(200 , 200) 
character* 10 FNOPEN 
FN0PEN=’ inp02.dat * 

OPEN (UNIT=1 , FILE=FNOPEN , F0RM= ’ FORMATTED ' , STATUS= ’ OLD ’ ) 
READ(1,*) JGRD.KGRD 

READ(1,*) ( (X( J ,K) , J-l , JGRD) ,K=1 ,KGRD) 

READ(1,*) ( (Y(J ,K) ,J=1,JGRD),K=1,KGRD) 

CLOSE (1) 



FN0PEN=’ inp02.dat’ 

OPEN (UNIT=2 , FILE=FNOPEN , F0RM= ’ FORMATTED ’ , STATUS= ’ UNKNOWN ’ ) 

WRITE(2, *) JGRD-1 ,KGRD 

WRITE (2 , * ) ( (X( J , K) , J=2 , JGRD) , K=1 ,KGRD ) 

WRITE(2, *) ((Y(J,K) , J=2, JGRD) ,K=1,KGRD) 

CL0SE(2) 



FNOPEN=’ inp03.dat’ 

OPEN (UNIT= 1 ,FILE=FNOPEN , FORM= ’ FORMATTED ’ , STATUS= ’ OLD ’ ) 
READ(1,*) JGRD ,KGRD 

READ(1,*) ((X(J,K),J=1,JGRD),K=1,KGRD) 

READ(1,*) ((Y(J,K),J=1, JGRD),K=1,KGRD) 

CLOSE (1) 



FNOPEN=’ inp03.dat’ 

0PEN(UNIT=2 , FILE=FNOPEN , FORM= ’ FORMATTED ’ , STATUS= ’ UNKNOWN ' ) 

WRITE(2 , *) JGRD-1 ,KGRD 

WRITE (2, *) ( (X(J ,K) ,J=2,JGRD) ,K=1,KGRD) 

WRITE(2, *) ((Y(J,K) > J=2,JGRD) > K=1,KGRD) 

CL0SE(2) 



FNOPEN=’ inp05.dat’ 

OPEN (UNIT=1 , FILE=FNOPEN , FORM= ' FORMATTED * , STATUS= ’ OLD ’ ) 
READ(1,*) JGRD ,KGRD 

READ ( 1 , * ) ((X(J,K),J=1,JGRD),K=1,KGRD) 

READ(1,*) ( (Y( J ,K) ,J=1,JGRD),K=1,KGRD) 

CLOSE(l) 

FNOPEN=’ inp05.dat’ 

OPEN (UNIT=2 , FILE=FNOPEN , FORM= ' FORMATTED ’ , STATUS^ ’ UNKNOWN ’ ) 

WRITE (2,*) JGRD-1 ,KGRD 

WRITE(2, *) ((X(J,K),J=2,JGRD),K=1,KGRD) 

WRITE(2 , *) ((Y(J,K),J=2,JGRD),K=1,KGRD) 

CL0SE(2) 



158 



FNOPEN= ’ inp06 . dat ’ 

0PEN(UNIT=1 , FI LE=FNO P EN , FORM= ’ FORMATTED ’ , STATUS= ’ OLD ’ ) 
READ(1,*) JGRD , KGRD 

READ ( 1 , * ) ((X(J,K),J=1,JGRD),K=1,KGRD) 

READ(1,*) ((Y(J, K),J=1, JGRD), K=l, KGRD) 

CLOSE(l) 



FNOPEN=' inp06.dat’ 

OPEN (UNIT=2 , FILE=FNOPEN, FORM= ’ FORMATTED ’ , STATUS= ’ UNKNOWN ' ) 

WRITE (2,*) JGRD- 1, KGRD 

WRITE(2,*) ((X(J,K) ,J=2, JGRD), K=l, KGRD) 

WRITE(2,*) ((Y(J,K),J=2, JGRD), K=l, KGRD) 

CL0SE(2) 



stop 

end 
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Appendix F 

Case Study Foil Offsets 



Offsets start at the trailing edge of the pressure side and march forward around the 
leading edge to the trailing edge of the suction side. Both data sets use a normalized 
chord length 1.0. The offsets listed here correspond to the foils displayed in Figures 
5-1 and 5-9. 
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F.l Case Study I: The HRA 


0.550 


0.047207 


Foil 


0.600 

0.650 

0.700 


0.045990 

0.044037 

0.041199 


1.0 0.0 


0.725 


0.039327 


0.995 - 0.000896 


0.750 


0.037072 


0.990 - 0.000797 


0.775 


0.034427 


0.985 - 0.000703 


0.800 


0.031406 


0.975 - 0.000537 


0.825 


0.028036 


0.965 - 0.000397 


0.850 


0.024371 


0.950 - 0.000225 


0.875 


0.020515 


0.925 - 0.000031 


0.900 


0.016571 


0.900 0.000049 


0.925 


0.012588 


0.875 0.000009 


0.950 


0.008599 


0.850 - 0.000150 


0.965 


0.006225 


0.825 - 0.000405 


0.975 


0.004671 


0.800 - 0.000734 


0.985 


0.003166 


0.775 - 0.001126 


0.990 


0.002433 


0.750 - 0.001574 


0.995 


0.001713 


0.725 - 0.002076 


1.000 


0.0 



0.700 - 0.002631 
0.650 - 0.003968 
0.600 - 0.005685 
0.550 - 0.007780 
0.500 - 0.010205 
0.450 - 0.012882 
0.400 - 0.015697 
0.350 - 0.018500 
0.300 - 0.021093 
0.250 - 0.023204 
0.200 - 0.024436 
0.175 - 0.024512 
0.150 - 0.024077 
0.125 - 0.023054 
0.100 - 0.021373 
0.075 - 0.018969 
0.050 - 0.015681 
0.035 - 0.013111 
0.025 - 0.011017 
0.015 - 0.008420 
0.010 - 0.006793 
0.005 - 0.004709 
0.0025 - 0.003274 
0.001 - 0.002036 
0.0 0.0 
0.001 0.002076 
0.0025 0.003375 
0.005 0.004911 
0.010 0.007200 
0.015 0.009036 
0.025 0.012061 
0.035 0.014599 
0.050 0.017865 
0.075 0.022406 
0.100 0.026204 
0.125 0.029474 
0.150 0.032329 
0.175 0.034841 
0.200 0.037059 
0.250 0.040736 
0.300 0.043546 
0.350 0.045599 
0.400 0.046962 
0.450 0.047677 
0.500 0.047760 
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2 B-l Foil With Cup Mod- 


0.069396 


- 0.012064 


ification To Trailing Edge 


0.060654 

0.053133 


- 0.011696 

- 0.011315 






0.046637 


- 0.010923 


1.000000 


0.000000 


0.040974 


- 0.010524 


0.989517 


0.002052 


0.035949 


- 0.010120 


0.979110 


0.003558 


0.031370 


- 0.009713 


0.968851 


0.004586 


0.027095 


- 0.009301 


0.958810 


0.005204 


0.023110 


- 0.008873 


0.949059 


0.005478 


0.019418 


- 0.008416 


0.939668 


0.005477 


0.016024 


- 0.007916 


0.930708 


0.005269 


0.012931 


- 0.007361 


0.922250 


0.004920 


0.010143 


- 0.006736 


0.914365 


0.004499 


0.007665 


- 0.006030 


0.907121 


0.004072 


0.005500 


- 0.005229 


0.900497 


0.003672 


0.003652 


- 0.004320 


0.894374 


0.003298 


0.002127 


- 0.003289 


0.888628 


0.002943 


0.000963 


- 0.002132 


0.883131 


0.002604 


0.000220 


- 0.000846 


0.877761 


0.002274 


0.000000 


0.000570 


0.872391 


0.001950 


0.000133 


0.001781 


0.866896 


0.001626 


0.000676 


0.003074 


0.861152 


0.001297 


0.001617 


0.004450 


0.855033 


0.000959 


0.002986 


0.005909 


0.848423 


0.000606 


0.004814 


0.007453 


0.841250 


0.000236 


0.007129 


0.009082 


0.833456 


- 0.000154 


0.009962 


0.010796 


0.824986 


- 0.000568 


0.013344 


0.012596 


0.815783 


- 0.001006 


0.017305 


0.014480 


0.805789 


- 0.001473 


0.021881 


0.016443 


0.794947 


- 0.001971 


0.027105 


0.018482 


0.783201 


- 0.002503 


0.033011 


0.020591 


0.770495 


- 0.003071 


0.039635 


0.022766 


0.756770 


- 0.003678 


0.047010 


0.025003 


0.742000 


- 0.004322 


0.055170 


0.027297 


0.726196 


- 0.004997 


0.064149 


0.029643 


0.709371 


- 0.005694 


0.073983 


0.032038 


0.691542 


- 0.006406 


0.084704 


0.034476 


0.672721 


- 0.007126 


0.096348 


0.036953 


0.652924 


- 0.007845 


0.108948 


0.039464 


0.632166 


- 0.008556 


0.122492 


0.042000 


0.610461 


- 0.009251 


0.136919 


0.044547 


0.587823 


- 0.009922 


0.152161 


0.047087 


0.564277 


- 0.010563 


0.168152 


0.049605 


0.539906 


- 0.011168 


0.184827 


0.052085 


0.514822 


- 0.011734 


0.202117 


0.054513 


0.489138 


- 0.012255 


0.219958 


0.056871 


0.462964 


- 0.012730 


0.238283 


0.059145 


0.436413 


- 0.013153 


0.257024 


0.061319 


0.409596 


- 0.013521 


0.276117 


0.063376 


0.382625 


- 0.013830 


0.295494 


0.065302 


0.355611 


- 0.014076 


0.315089 


0.067080 


0.328667 


- 0.014255 


0.334845 


0.068702 


0.301935 


- 0.014366 


0.354719 


0.070169 


0.275618 


- 0.014414 


0.374667 


0.071483 


0.249922 


- 0.014402 


0.394646 


0.072647 


0.225056 


- 0.014337 


0.414616 


0.073662 


0.201226 


- 0.014221 


0.434532 


0.074531 


0.178641 


- 0.014061 


0.454353 


0.075256 


0.157509 


- 0.013861 


0.474035 


0.075838 


0.138036 


- 0.013626 


0.493537 


0.076282 


0.120430 


- 0.013359 


0.512816 


0.076587 


0.104875 


- 0.013067 


0.531828 


0.076757 


0.091313 


- 0.012752 


0.550533 


0.076795 


0.079551 


- 0.012417 


0.568899 


0.076701 
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0.586921 

0.604596 

0.621920 

0.638890 

0.655504 

0.671758 

0.687648 

0.703172 

0.718327 

0.733108 

0.747514 

0.761542 

0.775192 

0.788463 

0.801358 

0.813876 

0.826017 

0.837782 

0.849171 

0.860186 

0.870825 

0.881089 

0.890979 

0.900510 

0.909739 

0.918732 

0.927555 

0.936276 

0.944959 

0.953672 

0.962481 

0.971452 

0.980651 

0.990145 

1.000000 



0.076477 

0.076125 

0.075647 

0.075044 

0.074317 

0.073469 

0.072500 

0.071413 

0.070208 

0.068888 

0.067454 

0.065910 

0.064268 

0.062543 

0.060746 

0.058891 

0.056991 

0.055059 

0.053108 

0.051152 

0.049203 

0.047275 

0.045380 

0.043511 

0.041590 

0.039530 

0.037241 

0.034634 

0.031619 

0.028109 

0.024014 

0.019244 

0.013711 

0.007326 

0.000000 
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